home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Source Code / Libraries / DCLAP 4j / SeqPups / apps / fastdnaml.src / fastDNAml_1_0_6.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-12-17  |  127.5 KB  |  4,682 lines  |  [TEXT/R*ch]
  1.  
  2. #if  defined(__MWERKS__) || defined(THINK_C)  || defined(MPW)
  3. #define MACINTOSH 1
  4. #endif
  5.  
  6. #define programName     "fastDNAml"
  7. #define programVersion  "1.0.6"
  8. #define programDate     "November 20, 1992"
  9.  
  10. /*  Copyright notice from dnaml:
  11.  *
  12.  *  version 3.3. (c) Copyright 1986, 1990 by the University of Washington
  13.  *  and Joseph Felsenstein.  Written by Joseph Felsenstein.  Permission is
  14.  *  granted to copy and use this program provided no fee is charged for it
  15.  *  and provided that this copyright notice is not removed.
  16.  *
  17.  *  Conversion to C and changes in sequential code by Gary Olsen, 1991-1992
  18.  *
  19.  *  p4 version by Hideo Matsuda and Ross Overbeek, 1991-1992
  20.  */
  21.  
  22. /*
  23.  *  1.0.1  Add ntaxa to tree comments
  24.  *         Set minimum branch length on reading tree
  25.  *         Add blanks around operators in treeString (for prolog parsing)
  26.  *         Add program version to treeString comments
  27.  *
  28.  *  1.0.2  Improved option line diagnostics
  29.  *         Improved auxiliary line diagnostics
  30.  *         Removed some trailing blanks from output
  31.  *
  32.  *  1.0.3  Checkpoint trees that do not need any optimization
  33.  *         Print restart tree likelihood before optimizing
  34.  *         Fix treefile option so that it really toggles
  35.  *
  36.  *  1.0.4  Add support for tree fact (instead of true Newick tree) in
  37.  *            processTreeCom, treeReadLen, str_processTreeCom and
  38.  *            str_treeReadLen
  39.  *         Use bit operations in randum
  40.  *         Correct error in bootstrap mask used with weighting mask
  41.  *
  42.  *  1.0.5  Fix reading of underscore as first nonblank character in name
  43.  *         Add strchr and strstr functions to source code
  44.  *         Add output treefile name to message "Tree also written ..."
  45.  *
  46.  *  1.0.6  Change (! nsites) test in setupTopol to (nsites == 0) for MIPS R4000
  47.  *         Add vectorizing compiler directives for CRAY
  48.  *         Include updates and corrections to parallel code from H. Matsuda
  49.  */
  50.  
  51. #ifdef Master
  52. #  undef  Master
  53. #  define Master     1
  54. #  define Slave      0
  55. #  define Sequential 0
  56. #else
  57. #  ifdef Slave
  58. #    undef Slave
  59. #    define Master     0
  60. #    define Slave      1
  61. #    define Sequential 0
  62. #  else
  63. #    ifdef Sequential
  64. #      undef Sequential
  65. #    endif
  66. #    define Master     0
  67. #    define Slave      0
  68. #    define Sequential 1
  69. #  endif
  70. #endif
  71.  
  72. #include <stdio.h>
  73. #include <math.h>
  74. #include "fastDNAml_1_0.h"
  75.  
  76. #if Master || Slave
  77. #include "p4.h"
  78. #include "comm_link.h" 
  79. #endif
  80.  
  81. /*  Global variables */
  82.  
  83. longer  seed;      /*  random number seed with 6 bits per element */
  84.  
  85. xarray
  86.     *usedxtip, *freextip;
  87.  
  88. static double
  89.     rate[maxcategories+1],    /*  rates for categories */
  90.     ratvalue[maxpatterns],    /*  rates per site */
  91.     wgt_rate[maxpatterns],    /*  weighted rate per site */
  92.     wgt_rate2[maxpatterns];   /*  weighted rate**2 per site */
  93.  
  94. static double
  95.     xi, xv, ttratio,          /*  transition/transversion info */
  96.     freqa, freqc, freqg, freqt,  /*  base frequencies */
  97.     freqr, freqy, invfreqr, invfreqy,
  98.     freqar, freqcy, freqgr, freqty,
  99.     totalwrate,    /*  total of weighted rate values */
  100.     fracchange;    /*  random matching fraquency (in a sense) */
  101.  
  102. static int
  103.     alias[maxsites+1],        /*  site corresponding to pattern */
  104.     aliasweight[maxsites+1],  /*  weight of given pattern */
  105.     category[maxsites+1],     /*  rate category of sequence positions */
  106.     catnumb[maxpatterns],     /*  category number by pattern number */
  107.     weight[maxsites+1];       /*  weight of sequence positions */
  108.  
  109. int
  110.     categs,        /*  number of rate categories */
  111.     endsite,       /*  number of unique sequence patterns */
  112.     extrainfo,     /*  runtime information switch */
  113.     numsp,         /*  number of species (same as tr->mxtips) */
  114.     nkeep,         /*  number of best trees to keep */
  115.     sites,         /*  number of input sequence positions */
  116.     trout,         /*  write tree to "treefile" */
  117.     weightsum;     /*  sum of weights of positions in analysis */
  118.  
  119. boolean
  120.     anerror,       /*  error flag */
  121.     bootstrap,     /*  do a set of bootstrap column weights */
  122.     freqsfrom,     /*  use empirical base frequencies */
  123.     interleaved,   /*  input data are in interleaved format */
  124.     jumble,        /*  use random addition order */
  125.     outgropt,      /*  use user-supplied outgroup */
  126.     printdata,     /*  echo data to output stream */
  127.     fastadd,       /*  test addition sites without global smoothing */
  128.     restart,       /*  continue sequential addition to partial tree */
  129.     treeprint;     /*  print tree to output stream */
  130.  
  131. char
  132.     *y[maxsp+1];   /*  sequence data array */
  133.  
  134. #if Sequential     /*  Use standard input */
  135. #  undef   DNAML_STEP
  136. #  define  DNAML_STEP  0
  137. #ifdef NoSTDIN
  138.   FILE *INFILE;
  139. #else
  140. #  define  INFILE  stdin
  141. #endif
  142. #endif
  143.  
  144. #if Master
  145. #  define MAX_SEND_AHEAD 400
  146.   char   *best_tr_recv = NULL;     /* these are used for flow control */
  147.   double  best_lk_recv;
  148.   int     send_ahead = 0;          /* number of outstanding sends */
  149.  
  150. #  ifdef DNAML_STEP
  151. #    define  DNAML_STEP  1
  152. #  endif
  153. #  define  INFILE  Seqf
  154. #  define  OUTFILE Outf
  155.   FILE *INFILE, *OUTFILE;
  156.   comm_block comm_slave;
  157. #endif
  158.  
  159. #if Slave
  160. #  undef   DNAML_STEP
  161. #  define  DNAML_STEP  0
  162. #  define  INFILE  Seqf
  163. #  define  OUTFILE Outf
  164.   FILE *INFILE, *OUTFILE;
  165.   comm_block comm_master;
  166. #endif
  167.  
  168. #if DNAML_STEP
  169.   int begin_step_time, end_step_time;
  170. #endif
  171.  
  172. #if Debug
  173.   FILE *debug;
  174. #endif
  175.  
  176. #ifdef  CRAY
  177. #  define  Vectorize
  178. #endif
  179.  
  180. #if DNAML_STEP
  181. #  define  REPORT_ADD_SPECS  p4_send(DNAML_ADD_SPECS, DNAML_HOST_ID, NULL, 0)
  182. #  define  REPORT_SEND_TREE  p4_send(DNAML_SEND_TREE, DNAML_HOST_ID, NULL, 0)
  183. #  define  REPORT_RECV_TREE  p4_send(DNAML_RECV_TREE, DNAML_HOST_ID, NULL, 0)
  184. #  define  REPORT_STEP_TIME \
  185.    {\
  186.        char send_buf[80]; \
  187.        end_step_time = p4_clock(); \
  188.        (void) sprintf(send_buf, "%d", end_step_time-begin_step_time); \
  189.        p4_send(DNAML_STEP_TIME, DNAML_HOST_ID, send_buf,strlen(send_buf)+1); \
  190.        begin_step_time = end_step_time; \
  191.    }
  192. #else
  193. #  define  REPORT_ADD_SPECS
  194. #  define  REPORT_SEND_TREE
  195. #  define  REPORT_RECV_TREE
  196. #  define  REPORT_STEP_TIME
  197. #endif
  198.  
  199. /*=======================================================================*/
  200. /*                              PROGRAM                                  */
  201. /*=======================================================================*/
  202. /*                    Best tree handling for dnaml                       */
  203. /*=======================================================================*/
  204.  
  205. /*  Tip value comparisons
  206.  *
  207.  *  Use void pointers to hide type from other routines.  Only tipValPtr and
  208.  *  cmpTipVal need to be changed to alter the nature of the values compared
  209.  *  (e.g., names instead of node numbers).
  210.  *
  211.  *    cmpTipVal(tipValPtr(nodeptr p), tipValPtr(nodeptr q)) == -1, 0 or 1.
  212.  *
  213.  *  This provides direct comparison of tip values (for example, for
  214.  *  definition of tr->start).
  215.  */
  216.  
  217.  
  218. void  *tipValPtr (p)  nodeptr  p; { return  (void *) & p->number; }
  219.  
  220.  
  221. int  cmpTipVal (v1, v2)
  222.     void  *v1, *v2;
  223.   { /* cmpTipVal */
  224.     int  i1, i2;
  225.  
  226.     i1 = *((int *) v1);
  227.     i2 = *((int *) v2);
  228.     return  (i1 < i2) ? -1 : ((i1 == i2) ? 0 : 1);
  229.   } /* cmpTipVal */
  230.  
  231.  
  232. /*  These are the only routines that need to UNDERSTAND topologies */
  233.  
  234.  
  235. topol  *setupTopol (maxtips, nsites)
  236.     int     maxtips, nsites;
  237.   { /* setupTopol */
  238.     topol   *tpl;        
  239.  
  240.     if (! (tpl = (topol *) malloc( (unsigned long)sizeof(topol))) || 
  241.         ! (tpl->links = (connptr) malloc((unsigned long) ((2 * maxtips - 3) *
  242.                                                      sizeof(connect)))) || 
  243.         (nsites && ! (tpl->log_f
  244.                 = (double *) malloc((unsigned long) (nsites * sizeof(double)))))) {
  245.       printf("ERROR: Unable to get topology memory");
  246.       tpl = (topol *) NULL;
  247.       }
  248.  
  249.     else {
  250.       if (nsites == 0)  tpl->log_f = (double *) NULL;
  251.       tpl->likelihood  = unlikely;
  252.       tpl->start       = (node *) NULL;
  253.       tpl->nextlink    = 0;
  254.       tpl->ntips       = 0;
  255.       tpl->nextnode    = 0;
  256.       tpl->opt_level   = 0;     /* degree of branch swapping explored */
  257.       tpl->scrNum      = 0;     /* position in sorted list of scores */
  258.       tpl->tplNum      = 0;     /* position in sorted list of trees */
  259.       tpl->log_f_valid = 0;     /* log_f value sites */
  260.       tpl->smoothed    = FALSE; /* branch optimization converged? */
  261.       }
  262.  
  263.     return  tpl;
  264.   } /* setupTopol */
  265.  
  266.  
  267. void  freeTopol (tpl)
  268.     topol  *tpl;
  269.   { /* freeTopol */
  270.     free((char *) tpl->links);
  271.     if (tpl->log_f)  free((char *) tpl->log_f);
  272.     free((char *) tpl);
  273.   } /* freeTopol */
  274.  
  275.  
  276. int  saveSubtree (p, tpl)
  277.     nodeptr  p;
  278.     topol   *tpl;
  279.     /*  Save a subtree in a standard order so that earlier branches
  280.      *  from a node contain lower value tips than do second branches from
  281.      *  the node.  This code works with arbitrary furcations in the tree.
  282.      */
  283.   { /* saveSubtree */
  284.     connptr  r, r0;
  285.     nodeptr  q, s;
  286.     int      t, t0, t1;
  287.  
  288.     r0 = tpl->links;
  289.     r = r0 + (tpl->nextlink)++;
  290.     r->p = p;
  291.     r->q = q = p->back;
  292.     r->z = p->z;
  293.     r->descend = 0;                     /* No children (yet) */
  294.  
  295.     if (q->tip) {
  296.       r->valptr = tipValPtr(q);         /* Assign value */
  297.       }
  298.  
  299.     else {                              /* Internal node, look at children */
  300.       s = q->next;                      /* First child */
  301.       do {
  302.         t = saveSubtree(s, tpl);        /* Generate child's subtree */
  303.  
  304.         t0 = 0;                         /* Merge child into list */
  305.         t1 = r->descend;
  306.         while (t1 && (cmpTipVal(r0[t1].valptr, r0[t].valptr) < 0)) {
  307.           t0 = t1;
  308.           t1 = r0[t1].sibling;
  309.           }
  310.         if (t0) r0[t0].sibling = t;  else  r->descend = t;
  311.         r0[t].sibling = t1;
  312.  
  313.         s = s->next;                    /* Next child */
  314.         } while (s != q);
  315.  
  316.       r->valptr = r0[r->descend].valptr;   /* Inherit first child's value */
  317.       }                                 /* End of internal node processing */
  318.  
  319.     return  r - r0;
  320.   } /* saveSubtree */
  321.  
  322.  
  323. nodeptr  minSubtreeTip (p0)
  324.     nodeptr  p0;
  325.   { /* minTreeTip */
  326.     nodeptr  minTip, p, testTip;
  327.  
  328.     if (p0->tip) return p0;
  329.  
  330.     p = p0->next;
  331.     minTip = minSubtreeTip(p->back);
  332.     while ((p = p->next) != p0) {
  333.       testTip = minSubtreeTip(p->back);
  334.       if (cmpTipVal(tipValPtr(testTip), tipValPtr(minTip)) < 0)
  335.         minTip = testTip;
  336.       }
  337.     return minTip;
  338.   } /* minTreeTip */
  339.  
  340.  
  341. nodeptr  minTreeTip (p)
  342.     nodeptr  p;
  343.   { /* minTreeTip */
  344.     nodeptr  minp, minpb;
  345.  
  346.     minp  = minSubtreeTip(p);
  347.     minpb = minSubtreeTip(p->back);
  348.     return cmpTipVal(tipValPtr(minp), tipValPtr(minpb)) < 0 ? minp : minpb;
  349.   } /* minTreeTip */
  350.  
  351.  
  352. void saveTree (tr, tpl)
  353.     tree   *tr;
  354.     topol  *tpl;
  355.     /*  Save a tree topology in a standard order so that first branches
  356.      *  from a node contain lower value tips than do second branches from
  357.      *  the node.  The root tip should have the lowest value of all.
  358.      */
  359.   { /* saveTree */
  360.     connptr  r;
  361.     double  *tr_log_f, *tpl_log_f;
  362.     int  i;
  363.  
  364.     tpl->nextlink = 0;                             /* Reset link pointer */
  365.     r = tpl->links + saveSubtree(minTreeTip(tr->start), tpl);  /* Save tree */
  366.     r->sibling = 0;
  367.  
  368.     tpl->likelihood = tr->likelihood;
  369.     tpl->start      = tr->start;
  370.     tpl->ntips      = tr->ntips;
  371.     tpl->nextnode   = tr->nextnode;
  372.     tpl->opt_level  = tr->opt_level;
  373.     tpl->smoothed   = tr->smoothed;
  374.  
  375.     if (tpl_log_f = tpl->log_f) {
  376.       tr_log_f  = tr->log_f;
  377.       i = tpl->log_f_valid = tr->log_f_valid;
  378.       while (--i >= 0)  *tpl_log_f++ = *tr_log_f++;
  379.       }
  380.     else {
  381.       tpl->log_f_valid = 0;
  382.       }
  383.   } /* saveTree */
  384.  
  385.  
  386. void restoreTree (tpl, tr)
  387.     topol  *tpl;
  388.     tree   *tr;
  389.   { /* restoreTree */
  390.     void  hookup();
  391.     void  initrav();
  392.  
  393.     connptr  r;
  394.     nodeptr  p, p0;
  395.     double  *tr_log_f, *tpl_log_f;
  396.     int  i;
  397.  
  398. /*  Clear existing connections */
  399.  
  400.     for (i = 1; i <= 2*(tr->mxtips) - 2; i++) {  /* Uses p = p->next at tip */
  401.       p0 = p = tr->nodep[i];
  402.       do {
  403.         p->back = (nodeptr) NULL;
  404.         p = p->next;
  405.         } while (p != p0);
  406.       }
  407.  
  408. /*  Copy connections from topology */
  409.  
  410.     for (r = tpl->links, i = 0; i < tpl->nextlink; r++, i++) {
  411.       hookup(r->p, r->q, r->z);
  412.       }
  413.  
  414.     tr->likelihood = tpl->likelihood;
  415.     tr->start      = tpl->start;
  416.     tr->ntips      = tpl->ntips;
  417.     tr->nextnode   = tpl->nextnode;
  418.     tr->opt_level  = tpl->opt_level;
  419.     tr->smoothed   = tpl->smoothed;
  420.  
  421.     if (tpl_log_f = tpl->log_f) {
  422.       tr_log_f = tr->log_f;
  423.       i = tr->log_f_valid = tpl->log_f_valid;
  424.       while (--i >= 0)  *tr_log_f++ = *tpl_log_f++;
  425.       }
  426.     else {
  427.       tr->log_f_valid = 0;
  428.       }
  429.  
  430.     initrav(tr->start);
  431.     initrav(tr->start->back);
  432.   } /* restoreTree */
  433.  
  434.  
  435. int initBestTree (bt, newKeep, nsites)
  436.     bestlist  *bt;
  437.     int        newKeep, nsites;
  438.   { /* initBestTree */
  439.     int  i, nlogf;
  440.  
  441.     if (newKeep < 1) newKeep = 1;
  442.     if (newKeep > maxkeep) newKeep = maxkeep;
  443.  
  444.     bt->best     = unlikely;
  445.     bt->worst    = unlikely;
  446.     bt->nvalid   = 0;
  447.     bt->ninit    = -1;
  448.     bt->numtrees = 0;
  449.     bt->improved = FALSE;
  450.     if (! (bt->start = setupTopol(numsp, nsites)))  return  0;
  451.  
  452.     for (i = bt->ninit + 1; i <= newKeep; i++) {
  453.       nlogf = (i <= maxlogf) ? nsites : 0;
  454.       if (! (bt->byScore[i] = setupTopol(numsp, nlogf)))  break;
  455.       bt->byTopol[i] = bt->byScore[i];
  456.       bt->ninit = i;
  457.       }
  458.  
  459.     return  (bt->nkeep = bt->ninit);
  460.   } /* initBestTree */
  461.  
  462.  
  463. int resetBestTree (bt)
  464.     bestlist  *bt;
  465.   { /* resetBestTree */
  466.     bt->best = unlikely;
  467.     bt->worst = unlikely;
  468.     bt->nvalid = 0;
  469.     bt->improved = FALSE;
  470.   } /* resetBestTree */
  471.  
  472.  
  473. void  freeBestTree(bt)
  474.     bestlist  *bt;
  475.   { /* freeBestTree */
  476.     while (bt->ninit >= 0)  freeTopol(bt->byScore[(bt->ninit)--]);
  477.     freeTopol(bt->start);
  478.   } /* freeBestTree */
  479.  
  480.  
  481. int  saveBestTree (bt, tr)
  482.     bestlist  *bt;
  483.     tree      *tr;
  484.   { /* saveBestTree */
  485.     topol  *tpl;
  486.     int  scrNum, tplNum;
  487.  
  488.     if ((bt->nvalid < bt->nkeep) || (tr->likelihood > bt->worst)) {
  489.       scrNum = 1;
  490.       tplNum = 1;
  491.       tpl = bt->byScore[scrNum];
  492.       saveTree(tr, tpl);
  493.       tpl->scrNum = scrNum;
  494.       tpl->tplNum = tplNum;
  495.  
  496.       bt->byTopol[tplNum] = tpl;
  497.       bt->nvalid          = 1;
  498.       bt->best            = tr->likelihood;
  499.       bt->worst           = tr->likelihood;
  500.       bt->improved        = TRUE;
  501.       }
  502.  
  503.     else
  504.       scrNum = 0;
  505.  
  506.     return  scrNum;
  507.   } /* saveBestTree */
  508.  
  509.  
  510. int  startOpt (bt, tr)
  511.     bestlist  *bt;
  512.     tree      *tr;
  513.   { /* startOpt */
  514.     int  scrNum;
  515.  
  516.     bt->nvalid = 0;
  517.     scrNum = saveBestTree (bt, tr);
  518.     bt->improved = FALSE;
  519.     return  scrNum;
  520.   } /* startOpt */
  521.  
  522.  
  523. int  setOptLevel (bt, level)
  524.     bestlist *bt;
  525.     int       level;
  526.   { /* setOptLevel */
  527.     int  scrNum;
  528.  
  529.     if (! bt->improved) {
  530.       bt->byScore[1]->opt_level = level;
  531.       scrNum = 1;
  532.       }
  533.     else
  534.       scrNum = 0;
  535.  
  536.     return  scrNum;
  537.   } /* setOptLevel */
  538.  
  539.  
  540. int  recallBestTree (bt, rank, tr)
  541.     bestlist  *bt;
  542.     int        rank;
  543.     tree      *tr;
  544.   { /* recallBestTree */
  545.     if (rank < 1)  rank = 1;
  546.     if (rank > bt->nvalid)  rank = bt->nvalid;
  547.     if (rank > 0)  restoreTree(bt->byScore[rank], tr);
  548.     return  rank;
  549.   } /* recallBestTree */
  550.  
  551. /*=======================================================================*/
  552. /*                       End of best tree routines                       */
  553. /*=======================================================================*/
  554.  
  555.  
  556. #if 0
  557.   void hang(msg) char *msg; {printf("Hanging around: %s\n", msg); while(1);}
  558. #endif
  559.  
  560.  
  561. void getnums (tr)
  562.     tree  *tr;
  563.     /* input number of species, number of sites */
  564.   { /* getnums */
  565.     printf("\n%s, version %s, %s\n\n",
  566.             programName,
  567.             programVersion,
  568.             programDate);
  569.     printf("Based on Joseph Felsenstein's\n\n");
  570.     printf("Nucleic acid sequence Maximum Likelihood method, ");
  571.     printf("version 3.3\n\n");
  572.  
  573.     if (fscanf(INFILE, "%d %d", &numsp, &sites) != 2) {
  574.       printf("ERROR: Problem reading number of species and sites\n");
  575.       anerror = TRUE;
  576.       return;
  577.       }
  578.     printf("%d Species, %d Sites\n\n", numsp, sites);
  579.  
  580.     if (numsp > maxsp) {
  581.       printf("TOO MANY SPECIES: adjust CONSTants\n");
  582.       anerror = TRUE;
  583.       }
  584.     else if (numsp < 4) {
  585.       printf("TOO FEW SPECIES\n");
  586.       anerror = TRUE;
  587.       }
  588.  
  589.     if (sites > maxsites) {
  590.       printf("TOO MANY SITES: adjust CONSTants\n");
  591.       anerror = TRUE;
  592.       }
  593.     else if (sites < 1) {
  594.       printf("TOO FEW SITES\n");
  595.       anerror = TRUE;
  596.       }
  597.  
  598.     tr->mxtips = numsp;
  599.   } /* getnums */
  600.  
  601.  
  602. boolean digit (ch) int ch; {return (ch >= '0' && ch <= '9'); }
  603.  
  604.  
  605. boolean white (ch) int ch; { return (ch == ' ' || ch == '\n' || ch == '\t'); }
  606.  
  607.  
  608. void uppercase (chptr)
  609.       int *chptr;
  610.     /* convert character to upper case -- either ASCII or EBCDIC */
  611.   { /* uppercase */
  612.     int  ch;
  613.  
  614.     ch = *chptr;
  615.     if ((ch >= 'a' && ch <= 'i') || (ch >= 'j' && ch <= 'r')
  616.                                  || (ch >= 's' && ch <= 'z'))
  617.       *chptr = ch + 'A' - 'a';
  618.   } /* uppercase */
  619.  
  620.  
  621. int base36 (ch)
  622.     int ch;
  623.   { /* base36 */
  624.     if      (ch >= '0' && ch <= '9') return (ch - '0');
  625.     else if (ch >= 'A' && ch <= 'I') return (ch - 'A' + 10);
  626.     else if (ch >= 'J' && ch <= 'R') return (ch - 'J' + 19);
  627.     else if (ch >= 'S' && ch <= 'Z') return (ch - 'S' + 28);
  628.     else if (ch >= 'a' && ch <= 'i') return (ch - 'a' + 10);
  629.     else if (ch >= 'j' && ch <= 'r') return (ch - 'j' + 19);
  630.     else if (ch >= 's' && ch <= 'z') return (ch - 's' + 28);
  631.     else return -1;
  632.   } /* base36 */
  633.  
  634.  
  635. int itobase36 (i)
  636.     int  i;
  637.   { /* itobase36 */
  638.     if      (i <  0) return '?';
  639.     else if (i < 10) return (i      + '0');
  640.     else if (i < 19) return (i - 10 + 'A');
  641.     else if (i < 28) return (i - 19 + 'J');
  642.     else if (i < 36) return (i - 28 + 'S');
  643.     else return '?';
  644.   } /* itobase36 */
  645.  
  646.  
  647. int findch (c)
  648.     int  c;
  649.   { /* findch */
  650.     int ch;
  651.  
  652.     while ((ch = getc(INFILE)) != EOF && ch != c) ;
  653.     return  ch;
  654.   } /* findch */
  655.  
  656.  
  657. #if Master || Slave
  658. int str_findch (treestrp, c)
  659.     char **treestrp;
  660.     int  c;
  661.   { /* str_findch */
  662.     int ch;
  663.  
  664.     while ((ch = *(*treestrp)++) != NULL && ch != c) ;
  665.     return  ch;
  666.   } /* str_findch */
  667. #endif
  668.  
  669.  
  670. void inputcategories ()
  671.     /* reads the categories for each site */
  672.   { /* inputcategories */
  673.     int  i, ch, ci;
  674.  
  675.     for (i = 1; i <= nmlngth; i++)  (void) getc(INFILE);
  676.     i = 1;
  677.     while (i <= sites) {
  678.       ch = getc(INFILE);
  679.       ci = base36(ch);
  680.       if (ci >= 0 && ci <= categs)
  681.         category[i++] = ci;
  682.       else if (! white(ch)) {
  683.         printf("ERROR: Bad category character (%c) at site %d\n", ch, i);
  684.         anerror = TRUE;
  685.         return;
  686.         }
  687.       }
  688.  
  689.     if (findch('\n') == EOF) {      /* skip to end of line */
  690.       printf("ERROR: Missing newline at end of category data\n");
  691.       anerror = TRUE;
  692.       }
  693.   } /* inputcategories */
  694.  
  695.  
  696. void inputweights ()
  697.     /* input the character weights 0, 1, 2 ... 9, A, B, ... Y, Z */
  698.   { /* inputweights */
  699.     int i, ch, wi;
  700.  
  701.     for (i = 2; i <= nmlngth; i++)  (void) getc(INFILE);
  702.     weightsum = 0;
  703.     i = 1;
  704.     while (i <= sites) {
  705.       ch = getc(INFILE);
  706.       wi = base36(ch);
  707.       if (wi >= 0)
  708.         weightsum += weight[i++] = wi;
  709.       else if (! white(ch)) {
  710.         printf("ERROR: Bad weight character: '%c'", ch);
  711.         printf("       Weights in dnaml must be a digit or a letter.\n");
  712.         anerror = TRUE;
  713.         return;
  714.         }
  715.       }
  716.  
  717.     if (findch('\n') == EOF) {      /* skip to end of line */
  718.       printf("ERROR: Missing newline at end of weight data\n");
  719.       anerror = TRUE;
  720.       }
  721.   } /* inputweights */
  722.  
  723.  
  724. void getoptions (tr)
  725.     tree  *tr;
  726.   { /* getoptions */
  727.     longer  bseed;
  728.     long    inseed, binseed;    /*  random number seed input from user */
  729.     int     ch, i, j, extranum;
  730.     double  randum();
  731.  
  732.     bootstrap    =  FALSE;  /* Don't bootstrap column weights */
  733.     categs       =      0;  /* Number of rate categories */
  734.     extrainfo    =      0;  /* No extra runtime info unless requested */
  735.     freqsfrom    =  FALSE;  /* Use empirical base frequencies */
  736.     tr->global   =     -1;  /* Default search locale for optimum */
  737.     tr->partswap =      1;  /* Default to swap locally after insert */
  738.     interleaved  =   TRUE;  /* By default, data format is interleaved */
  739.     jumble       =  FALSE;  /* Use random addition sequence */
  740.     nkeep        =      1;  /* Keep only the one best tree */
  741.     tr->outgr    =      1;  /* Outgroup number */
  742.     outgropt     =  FALSE;  /* User-defined outgroup rooting */
  743.     printdata    =  FALSE;  /* Don't echo data to output stream */
  744.     fastadd      = Master;  /* Smooth branches globally in add */
  745.     rate[1]      =    1.0;  /* Rate values */
  746.     restart      =  FALSE;  /* Restart from user tree */
  747.     treeprint    =   TRUE;  /* Print tree to output stream */
  748.     trout        =      0;  /* Output tree file */
  749.     ttratio      =    2.0;  /* Transition/transversion rate ratio */
  750.     tr->userlen  =  FALSE;  /* User-supplied branch lengths */
  751.     tr->usertree =  FALSE;  /* User-defined tree topologies */
  752.     tr->userwgt  =  FALSE;  /* User-defined position weights */
  753.     extranum     =      0;
  754.  
  755.     while ((ch = getc(INFILE)) != '\n' && ch != EOF) {
  756.       uppercase(& ch);
  757.       switch (ch) {
  758.           case '1' : printdata    = ! printdata; break;
  759.           case '3' : treeprint    = ! treeprint; break;
  760.           case '4' : trout        = 1; break;
  761.           case 'B' : bootstrap    = TRUE; binseed = 0; extranum++; break;
  762.           case 'C' : categs       = -1; extranum++; break;
  763.           case 'E' : extrainfo    = -1; break;
  764.           case 'F' : freqsfrom    = TRUE; break;
  765.           case 'G' : tr->global   = -2; break;
  766.           case 'I' : interleaved  = ! interleaved; break;
  767.           case 'J' : jumble       = TRUE; inseed = 0; extranum++; break;
  768.           case 'K' : extranum++;          break;
  769.           case 'L' : tr->userlen  = TRUE; break;
  770.           case 'O' : outgropt     = TRUE; tr->outgr = 0; extranum++; break;
  771.           case 'Q' : fastadd      = ! Master; break;
  772.           case 'R' : restart      = TRUE; break;
  773.           case 'T' : ttratio      = -1.0; extranum++; break;
  774.           case 'U' : tr->usertree = TRUE; break;
  775.           case 'W' : tr->userwgt  = TRUE; weightsum = 0; extranum++; break;
  776.           case 'Y' : trout        = 1 - trout; break;
  777.           case ' ' : break;
  778.           case '\t': break;
  779.           default  :
  780.               printf("ERROR: Bad option character: '%c'\n", ch);
  781.               anerror = TRUE;
  782.               return;
  783.           }
  784.       }
  785.  
  786.     if (ch == EOF) {
  787.       printf("ERROR: End-of-file in options list\n");
  788.       anerror = TRUE;
  789.       return;
  790.       }
  791.  
  792.     if (tr->usertree && restart) {
  793.       printf("ERROR:  The restart and user-tree options conflict:\n");
  794.       printf("        Restart adds rest of taxa to a starting tree;\n");
  795.       printf("        User-tree does not add any taxa.\n\n");
  796.       anerror = TRUE;
  797.       return;
  798.       }
  799.     if (tr->usertree && jumble) {
  800.       printf("WARNING:  The jumble and user-tree options conflict:\n");
  801.       printf("          Jumble adds taxa to a tree in random order;\n");
  802.       printf("          User-tree does not use taxa addition.\n");
  803.       printf("          Jumble ingnored for this run.\n\n");
  804.       jumble = FALSE;
  805.       }
  806.     if (tr->userlen && tr->global != -1) {
  807.       printf("ERROR:  The global and user-lengths options conflict:\n");
  808.       printf("        Global optimizes a starting tree;\n");
  809.       printf("        User-lengths constrain the starting tree.\n\n");
  810.       anerror = TRUE;
  811.       return;
  812.       }
  813.     if (tr->userlen && ! tr->usertree) {
  814.       printf("WARNING:  User lengths required user tree option.\n");
  815.       printf("          User-tree option set for this run.\n\n");
  816.       tr->usertree = TRUE;
  817.       }
  818.  
  819.     /*  process lines with auxiliary data */
  820.  
  821.     while (extranum--) {
  822.       ch = getc(INFILE);
  823.       uppercase(& ch);
  824.       switch (ch) {
  825.  
  826.         case 'B':   /*  Bootstrap  */
  827.             if (! bootstrap || binseed != 0) {
  828.               printf("ERROR: Unexpected Bootstrap auxiliary data line\n");
  829.               anerror = TRUE;
  830.               }
  831.             else if (fscanf(INFILE,"%ld",&binseed)!=1 || findch('\n')==EOF) {
  832.               printf("ERROR: Problem reading boostrap random seed value\n");
  833.               anerror = TRUE;
  834.               }
  835.             break;
  836.  
  837.         case 'C':   /*  Categories  */
  838.             if (categs >= 0) {
  839.               printf("ERROR: Unexpected Categories auxiliary data line\n");
  840.               anerror = TRUE;
  841.               }
  842.             else if (fscanf(INFILE, "%d", &categs) != 1) {
  843.               printf("ERROR: Problem reading number of rate categories\n");
  844.               anerror = TRUE;
  845.               }
  846.             else if (categs < 1 || categs > maxcategories) {
  847.               printf("ERROR: Bad number of categories: %d\n", categs);
  848.               printf("Must be in range 1 - %d\n", maxcategories);
  849.               anerror = TRUE;
  850.               }
  851.             else {
  852.               for (j = 1; j <= categs
  853.                        && fscanf(INFILE, "%lf", &(rate[j])) == 1; j++) ;
  854.  
  855.               if ((j <= categs) || (findch('\n') == EOF)) {
  856.                 printf("ERROR: Problem reading rate values\n");
  857.                 anerror = TRUE;
  858.                 }
  859.               else
  860.                 inputcategories();
  861.               }
  862.             break;
  863.  
  864.         case 'E':   /*  Extra output information */
  865.             if (fscanf(INFILE,"%d",&extrainfo) != 1 || findch('\n') == EOF) {
  866.               printf("ERROR: Problem reading extra info value\n");
  867.               anerror = TRUE;
  868.               }
  869.             extranum++;  /*  Don't count this line since it is optional */
  870.             break;
  871.  
  872.         case 'G':   /*  Global  */
  873.             if (tr->global != -2) {
  874.               printf("ERROR: Unexpected Global auxiliary data line\n");
  875.               anerror = TRUE;
  876.               }
  877.             else if (fscanf(INFILE, "%d", &(tr->global)) != 1) {
  878.               printf("ERROR: Problem reading rearrangement region size\n");
  879.               anerror = TRUE;
  880.               }
  881.             else if (tr->global < 0) {
  882.               printf("WARNING: Global region size too small;\n");
  883.               printf("         value reset to local\n\n");
  884.               tr->global = 1;
  885.               }
  886.             else if (tr->global == 0)  tr->partswap = 0;
  887.             else if (tr->global > tr->mxtips - 3) {
  888.               tr->global = tr->mxtips - 3;
  889.               }
  890.  
  891.             while ((ch = getc(INFILE)) != '\n') {  /* Scan for second value */
  892.               if (! white(ch)) {
  893.                 if (ch != EOF)  (void) ungetc(ch, INFILE);
  894.                 if (ch == EOF || fscanf(INFILE, "%d", &(tr->partswap)) != 1
  895.                               || findch('\n') == EOF) {
  896.                   printf("ERROR: Problem reading insert swap region size\n");
  897.                   anerror = TRUE;
  898.                   }
  899.                 else if (tr->partswap < 0)  tr->partswap = 1;
  900.                 else if (tr->partswap > tr->mxtips - 3) {
  901.                   tr->partswap = tr->mxtips - 3;
  902.                   }
  903.  
  904.                 if (tr->partswap > tr->global)  tr->global = tr->partswap;
  905.                 break;   /*  Break while loop */
  906.                 }
  907.               }
  908.  
  909.             extranum++;  /*  Don't count this line since it is optional */
  910.             break;       /*  Break switch statement */
  911.  
  912.         case 'J':   /*  Jumble  */
  913.             if (! jumble || inseed != 0) {
  914.               printf("ERROR: Unexpected Jumble auxiliary data line\n");
  915.               anerror = TRUE;
  916.               }
  917.             else if (fscanf(INFILE,"%ld",&inseed) != 1 || findch('\n')==EOF) {
  918.               printf("ERROR: Problem reading jumble random seed value\n");
  919.               anerror = TRUE;
  920.               }
  921.             break;
  922.  
  923.         case 'K':   /*  Keep  */
  924.             if (fscanf(INFILE, "%d", &nkeep) != 1 || findch('\n') == EOF ||
  925.                 nkeep < 1) {
  926.               printf("ERROR: Problem reading number of kept trees\n");
  927.               anerror = TRUE;
  928.               }
  929.             else if (nkeep > maxkeep) {
  930.               printf("WARNING: Kept trees lowered from %d to %d\n\n",
  931.                       nkeep, maxkeep);
  932.               nkeep = maxkeep;
  933.               }
  934.             break;
  935.  
  936.         case 'O':   /*  Outgroup  */
  937.             if (! outgropt || tr->outgr > 0) {
  938.               printf("ERROR: Unexpected Outgroup auxiliary data line\n");
  939.               anerror = TRUE;
  940.               }
  941.             else if (fscanf(INFILE, "%d", &(tr->outgr)) != 1 ||
  942.                 findch('\n') == EOF) {
  943.               printf("ERROR: Problem reading outgroup number\n");
  944.               anerror = TRUE;
  945.               }
  946.             else if ((tr->outgr < 1) || (tr->outgr > tr->mxtips)) {
  947.               printf("ERROR: Bad outgroup: '%d'\n", tr->outgr);
  948.               anerror = TRUE;
  949.               }
  950.             break;
  951.  
  952.         case 'T':   /*  Transition/transversion ratio  */
  953.             if (ttratio >= 0.0) {
  954.               printf("ERROR: Unexpected Transition/transversion auxiliary data\n");
  955.               anerror = TRUE;
  956.               }
  957.             else if (fscanf(INFILE,"%lf",&ttratio)!=1 || findch('\n')==EOF) {
  958.               printf("ERROR: Problem reading transition/transversion ratio\n");
  959.               anerror = TRUE;
  960.               }
  961.             break;
  962.  
  963.         case 'W':    /*  Weights  */
  964.             if (! tr->userwgt || weightsum > 0) {
  965.               printf("ERROR: Unexpected Weights auxiliary data\n");
  966.               anerror = TRUE;
  967.               }
  968.             else {
  969.               inputweights();
  970.               }
  971.             break;
  972.  
  973.         case 'Y':    /*  Output tree file  */
  974.             if (! trout) {
  975.               printf("ERROR: Unexpected Treefile auxiliary data\n");
  976.               anerror = TRUE;
  977.               }
  978.             else if (fscanf(INFILE,"%d",&trout) != 1 || findch('\n') == EOF) {
  979.               printf("ERROR: Problem reading output tree-type number\n");
  980.               anerror = TRUE;
  981.               }
  982.             else if ((trout < 0) || (trout > 2)) {
  983.               printf("ERROR: Bad output tree-type number: '%d'\n", trout);
  984.               anerror = TRUE;
  985.               }
  986.             extranum++;  /*  Don't count this line since it is optional */
  987.             break;
  988.  
  989.         default:
  990.             printf("ERROR: Auxiliary options line starts with '%c'\n", ch);
  991.             anerror = TRUE;
  992.             break;
  993.         }
  994.  
  995.         if (anerror)  return;
  996.       }
  997.  
  998.     if (anerror) return;
  999.  
  1000.     if (! tr->userwgt) {
  1001.       for (i = 1; i <= sites; i++) weight[i] = 1;
  1002.       weightsum = sites;
  1003.       }
  1004.  
  1005.     if (tr->userwgt && weightsum < 1) {
  1006.       printf("ERROR:  Missing or bad user-supplied weight data.\n");
  1007.       anerror = TRUE;
  1008.       return;
  1009.       }
  1010.  
  1011.     if (bootstrap) {
  1012.       int  nonzero;
  1013.  
  1014.       printf("Bootstrap random number seed = %ld\n\n", binseed);
  1015.  
  1016.       if (binseed < 0)  binseed = -binseed;
  1017.       for (i = 0; i <= 5; i++) {bseed[i] = binseed & 63; binseed >>= 6;}
  1018.  
  1019.       nonzero = 0;
  1020.       for (i = 1; i <= sites; i++)  if (weight[i] > 0) nonzero++;
  1021.  
  1022.       for (j = 1; j <= nonzero; j++) aliasweight[j] = 0;
  1023.       for (j = 1; j <= nonzero; j++)
  1024.         aliasweight[(int) (nonzero*randum(bseed)) + 1]++;
  1025.  
  1026.       j = 0;
  1027.       weightsum = 0;
  1028.       for (i = 1; i <= sites; i++) {
  1029.         if (weight[i] > 0) {
  1030.           weightsum += (weight[i] *= aliasweight[++j]);
  1031.           }
  1032.         }
  1033.       }
  1034.  
  1035.     if (jumble) {
  1036.       printf("Jumble random number seed = %ld\n\n", inseed);
  1037.       if (inseed < 0)  inseed = -inseed;
  1038.       for (i = 0; i <= 5; i++) {seed[i] = inseed & 63; inseed >>= 6;}
  1039.       }
  1040.  
  1041.     if (categs > 0) {
  1042.       printf("Site category   Rate of change\n\n");
  1043.       for (i = 1; i <= categs; i++)
  1044.         printf("           %c%13.3f\n", itobase36(i), rate[i]);
  1045.       putchar('\n');
  1046.       for (i = 1; i <= sites; i++) {
  1047.         if ((weight[i] > 0) && (category[i] < 1)) {
  1048.           printf("ERROR: Bad category (%c) at site %d\n",
  1049.                   itobase36(category[i]), i);
  1050.           anerror = TRUE;
  1051.           return;
  1052.           }
  1053.         }
  1054.       }
  1055.     else if (categs < 0) {
  1056.       printf("ERROR: Category auxiliary data missing from input\n");
  1057.       anerror = TRUE;
  1058.       return;
  1059.       }
  1060.     else {                                              /* categs == 0 */
  1061.       for (i = 1; i <= sites; i++) category[i] = 1;
  1062.       categs = 1;
  1063.       }
  1064.  
  1065.     if (tr->outgr < 1) {
  1066.       printf("ERROR: Outgroup auxiliary data missing from input\n");
  1067.       anerror = TRUE;
  1068.       return;
  1069.       }
  1070.  
  1071.     if (ttratio < 0.0) {
  1072.       printf("ERROR: Transition/transversion auxiliary data missing from input\n");
  1073.       anerror = TRUE;
  1074.       return;
  1075.       }
  1076.  
  1077.     if (tr->global < 0) {
  1078.       if (tr->global == -2) tr->global = tr->mxtips - 3;  /* Default global */
  1079.       else                  tr->global = tr->usertree ? 0 : 1; /* No global */
  1080.       }
  1081.  
  1082.     if (restart) {
  1083.       printf("Restart option in effect.  ");
  1084.       printf("Sequence addition will start from appended tree.\n\n");
  1085.       }
  1086.  
  1087.     if (tr->usertree && ! tr->global) {
  1088.       printf("User-supplied tree topology%swill be used.\n\n",
  1089.         tr->userlen ? " and branch lengths " : " ");
  1090.       }
  1091.     else {
  1092.       if (! tr->usertree) {
  1093.         printf("Rearrangements of partial trees may cross %d %s.\n",
  1094.                tr->partswap, tr->partswap == 1 ? "branch" : "branches");
  1095.         }
  1096.       printf("Rearrangements of full tree may cross %d %s.\n\n",
  1097.              tr->global, tr->global == 1 ? "branch" : "branches");
  1098.       }
  1099.   } /* getoptions */
  1100.  
  1101.  
  1102. void getbasefreqs ()
  1103.   { /* getbasefreqs */
  1104.     double  suma, sumb;
  1105.  
  1106.     if (freqsfrom) printf("Empirical ");
  1107.     printf("Base Frequencies:\n\n");
  1108.  
  1109.     if (! freqsfrom) {
  1110.       if (fscanf(INFILE, "%lf%lf%lf%lf", &freqa, &freqc, &freqg, &freqt) != 4
  1111.                                                 || findch('\n') == EOF) {
  1112.         printf("ERROR: Problem reading user base frequencies\n");
  1113.         anerror = TRUE;
  1114.         return;
  1115.         }
  1116.       }
  1117.  
  1118.     printf("   A    %10.5f\n", freqa);
  1119.     printf("   C    %10.5f\n", freqc);
  1120.     printf("   G    %10.5f\n", freqg);
  1121.     printf("  T(U)  %10.5f\n\n", freqt);
  1122.     freqr = freqa + freqg;
  1123.     invfreqr = 1.0/freqr;
  1124.     freqar = freqa * invfreqr;
  1125.     freqgr = freqg * invfreqr;
  1126.     freqy = freqc + freqt;
  1127.     invfreqy = 1.0/freqy;
  1128.     freqcy = freqc * invfreqy;
  1129.     freqty = freqt * invfreqy;
  1130.     printf("Transition/transversion ratio = %10.6f\n\n", ttratio);
  1131.     suma = ttratio*freqr*freqy - (freqa*freqg + freqc*freqt);
  1132.     sumb = freqa*freqgr + freqc*freqty;
  1133.     xi = suma/(suma+sumb);
  1134.     xv = 1.0 - xi;
  1135.     if (xi <= 0.0) {
  1136.       printf("WARNING: This transition/transversion ratio\n");
  1137.       printf("         is impossible with these base frequencies!\n");
  1138.       printf("Transition/transversion parameter reset\n\n");
  1139.       xi = 3/5;
  1140.       xv = 2/5;
  1141.       }
  1142.     printf("(Transition/transversion parameter = %10.6f)\n\n", xi/xv);
  1143.     fracchange = 2.0 * xi*(freqa*freqgr + freqc*freqty)
  1144.                      + xv*(1.0 - freqa*freqa - freqc*freqc
  1145.                                - freqg*freqg - freqt*freqt);
  1146.     /*  fracchange = 2.0 * (xi * (freqa*freqgr + freqc*freqty)
  1147.                             xv * (freqa*freqg + freqc*freqt + freqr*freqy)
  1148.                            )
  1149.                    = 2.0 * (xi * sumb
  1150.                             xv * (freqa*freqg + freqc*freqt + freqr*freqy)
  1151.                            )
  1152.      */
  1153.   } /* getbasefreqs */
  1154.  
  1155.  
  1156. void getyspace ()
  1157.   { /* getyspace */
  1158.     long size;
  1159.     int  i;
  1160.     char *y0;
  1161.         unsigned long ysize;
  1162.         
  1163.     size = 4 * (sites/4 + 1);
  1164.     ysize= (numsp+1) * size * sizeof(char);
  1165.       y0 = malloc(ysize);
  1166.     if (! y0) {
  1167.       printf("ERROR: Unable to obtain space for data array\n");
  1168.       anerror = TRUE;
  1169.       return;
  1170.       }
  1171.  
  1172.     for (i = 0; i <= numsp; i++) {
  1173.       y[i] = y0;
  1174.       y0 += size;
  1175.       }
  1176.   } /* getyspace */
  1177.  
  1178.  
  1179. void setupTree (tr)
  1180.     tree  *tr;
  1181.   { /* setupTree */
  1182.     nodeptr  p0, p, q;
  1183.     int  i, j, tips, inter;
  1184.         unsigned long msize;
  1185.         
  1186.     tips  = tr->mxtips;
  1187.     inter = tr->mxtips - 1;
  1188.  
  1189.     msize= (tips + 3*inter) * sizeof(node);
  1190.       p0 = (nodeptr) malloc(msize);
  1191.     if (!p0) {
  1192.       printf("ERROR: Unable to obtain sufficient tree memory\n");
  1193.       anerror = TRUE;
  1194.       return;
  1195.       }
  1196.  
  1197.     tr->nodep[0] = (node *) NULL;    /* Use as 1-based array */
  1198.  
  1199.     for (i = 1; i <= tips; i++) {   /* Set-up tips */
  1200.       p = p0++;
  1201.       p->x      = (xarray *) NULL;
  1202.       p->tip    = (char *) NULL;
  1203.       p->number = i;
  1204.       p->next   = p;
  1205.       p->back   = (node *) NULL;
  1206.  
  1207.       tr->nodep[i] = p;
  1208.       }
  1209.  
  1210.     for (i = tips + 1; i <= tips + inter; i++) { /* Internal nodes */
  1211.       q = (node *) NULL;
  1212.       for (j = 1; j <= 3; j++) {
  1213.         p = p0++;
  1214.         p->x      = (xarray *) NULL;
  1215.         p->tip    = (char *) NULL;
  1216.         p->number = i;
  1217.         p->next   = q;
  1218.         p->back   = (node *) NULL;
  1219.         q = p;
  1220.         }
  1221.       p->next->next->next = p;
  1222.       tr->nodep[i] = p;
  1223.       }
  1224.  
  1225.     tr->likelihood  = unlikely;
  1226.     tr->start       = (node *) NULL;
  1227.     tr->outgrnode   = tr->nodep[tr->outgr];
  1228.     tr->ntips       = 0;
  1229.     tr->nextnode    = 0;
  1230.     tr->opt_level   = 0;
  1231.     tr->smoothed    = FALSE;
  1232.     tr->log_f_valid = 0;
  1233.   } /* setupTree */
  1234.  
  1235.  
  1236. void freeTreeNode (p)       /* Free tree node (sector) associated data */
  1237.     nodeptr  p;
  1238.   { /* freeTreeNode */
  1239.     if (p) {
  1240.       if (p->x) {
  1241.         if (p->x->a) free((char *) p->x->a);
  1242.         free((char *) p->x);
  1243.         }
  1244.       }
  1245.   } /* freeTreeNode */
  1246.  
  1247.  
  1248. void freeTree (tr)
  1249.     tree  *tr;
  1250.   { /* freeTree */
  1251.     nodeptr  p, q;
  1252.     int  i, tips, inter;
  1253.  
  1254.     tips  = tr->mxtips;
  1255.     inter = tr->mxtips - 1;
  1256.  
  1257.     for (i = 1; i <= tips; i++) freeTreeNode(tr->nodep[i]);
  1258.  
  1259.     for (i = tips + 1; i <= tips + inter; i++) {
  1260.       if (p = tr->nodep[i]) {
  1261.         if (q = p->next) {
  1262.           freeTreeNode(q->next);
  1263.           freeTreeNode(q);
  1264.           }
  1265.         freeTreeNode(p);
  1266.         }
  1267.       }
  1268.  
  1269.     free((char *) tr->nodep[1]);       /* Free the actual nodes */
  1270.   } /* freeTree */
  1271.  
  1272.  
  1273. void getdata (tr)
  1274.     tree  *tr;
  1275.     /* read sequences */
  1276.   { /* getdata */
  1277.     int  i, j, k, l, basesread, basesnew, ch;
  1278.     int  meaning[256];          /*  meaning of input characters */ 
  1279.     char *nameptr;
  1280.     boolean  allread, firstpass;
  1281.  
  1282.     for (i = 0; i <= 255; i++) meaning[i] = 0;
  1283.     meaning['A'] =  1;
  1284.     meaning['B'] = 14;
  1285.     meaning['C'] =  2;
  1286.     meaning['D'] = 13;
  1287.     meaning['G'] =  4;
  1288.     meaning['H'] = 11;
  1289.     meaning['K'] = 12;
  1290.     meaning['M'] =  3;
  1291.     meaning['N'] = 15;
  1292.     meaning['O'] = 15;
  1293.     meaning['R'] =  5;
  1294.     meaning['S'] =  6;
  1295.     meaning['T'] =  8;
  1296.     meaning['U'] =  8;
  1297.     meaning['V'] =  7;
  1298.     meaning['W'] =  9;
  1299.     meaning['X'] = 15;
  1300.     meaning['Y'] = 10;
  1301.     meaning['?'] = 15;
  1302.     meaning['-'] = 15;
  1303.  
  1304.     basesread = basesnew = 0;
  1305.  
  1306.     allread = FALSE;
  1307.     firstpass = TRUE;
  1308.     ch = ' ';
  1309.  
  1310.     while (! allread) {
  1311.       for (i = 1; i <= tr->mxtips; i++) {     /*  Read data line */
  1312.  
  1313.         if (firstpass) {                      /*  Read species names */
  1314.           j = 1;
  1315.           while (white(ch = getc(INFILE))) {  /*  Skip blank lines */
  1316.             if (ch == '\n')  j = 1;  else  j++;
  1317.             }
  1318.  
  1319.           if (j > nmlngth) {
  1320.             printf("ERROR: Blank name for species %d; ", i);
  1321.             printf("check number of species,\n");
  1322.             printf("       number of sites, and interleave option.\n");
  1323.             anerror = TRUE;
  1324.             return;
  1325.             }
  1326.  
  1327.           nameptr = tr->nodep[i]->name;
  1328.           for (k = 1; k < j; k++)  *nameptr++ = ' ';
  1329.  
  1330.           while (ch != '\n' && ch != EOF) {
  1331.             if (ch == '_' || white(ch))  ch = ' ';
  1332.             *nameptr++ = ch;
  1333.             if (++j > nmlngth) break;
  1334.             ch = getc(INFILE);
  1335.             }
  1336.  
  1337.           while (j++ <= nmlngth) *nameptr++ = ' ';
  1338.           *nameptr = '\0';                      /*  add null termination */
  1339.  
  1340.           if (ch == EOF) {
  1341.             printf("ERROR: End-of-file in name of species %d\n", i);
  1342.             anerror = TRUE;
  1343.             return;
  1344.             }
  1345.           }    /* if (firstpass) */
  1346.  
  1347.         j = basesread;
  1348.         while ((j < sites) && ((ch = getc(INFILE)) != EOF)
  1349.                            && ((! interleaved) || (ch != '\n'))) {
  1350.           uppercase(& ch);
  1351.           if (meaning[ch] || ch == '.') {
  1352.             j++;
  1353.             if (ch == '.') {
  1354.               if (i != 1) ch = y[1][j];
  1355.               else {
  1356.                 printf("ERROR: Dot (.) found at site %d of sequence 1\n", j);
  1357.                 anerror = TRUE;
  1358.                 return;
  1359.                 }
  1360.               }
  1361.             y[i][j] = ch;
  1362.             }
  1363.           else if (white(ch) || digit(ch)) ;
  1364.           else {
  1365.             printf("ERROR: Bad base (%c) at site %d of sequence %d\n",
  1366.                     ch, j, i);
  1367.             anerror = TRUE;
  1368.             return;
  1369.             }
  1370.           }
  1371.  
  1372.         if (ch == EOF) {
  1373.           printf("ERROR: End-of-file at site %d of sequence %d\n", j, i);
  1374.           anerror = TRUE;
  1375.           return;
  1376.           }
  1377.  
  1378.         if (! firstpass && (j == basesread)) i--;        /* no data on line */
  1379.         else if (i == 1) basesnew = j;
  1380.         else if (j != basesnew) {
  1381.           printf("ERROR: Sequences out of alignment\n");
  1382.           printf("%d (instead of %d) residues read in sequence %d\n",
  1383.                   j - basesread, basesnew - basesread, i);
  1384.           anerror = TRUE;
  1385.           return;
  1386.           }
  1387.  
  1388.         while (ch != '\n' && ch != EOF) ch = getc(INFILE);  /* flush line */
  1389.         }                                                  /* next sequence */
  1390.       firstpass = FALSE;
  1391.       basesread = basesnew;
  1392.       allread = (basesread >= sites);
  1393.       }
  1394.  
  1395.     /*  Print listing of sequence alignment */
  1396.  
  1397.     if (printdata) {
  1398.       j = nmlngth - 5 + ((sites + ((sites-1)/10))/2);
  1399.       if (j < nmlngth - 1) j = nmlngth - 1;
  1400.       if (j > 37) j = 37;
  1401.       printf("Name"); for (i=1;i<=j;i++) putchar(' '); printf("Sequences\n");
  1402.       printf("----"); for (i=1;i<=j;i++) putchar(' '); printf("---------\n");
  1403.       putchar('\n');
  1404.  
  1405.       for (i = 1; i <= sites; i += 60) {
  1406.         l = i + 59;
  1407.         if (l > sites) l = sites;
  1408.  
  1409.         if (tr->userwgt) {
  1410.           printf("Weights   ");
  1411.           for (j = 11; j <= nmlngth+3; j++) putchar(' ');
  1412.           for (k = i; k <= l; k++) {
  1413.             putchar(itobase36(weight[k]));
  1414.             if (((k % 10) == 0) && (k < l)) putchar(' ');
  1415.             }
  1416.           putchar('\n');
  1417.           }
  1418.  
  1419.         if (categs > 1) {
  1420.           printf("Categories");
  1421.           for (j = 11; j <= nmlngth+3; j++) putchar(' ');
  1422.           for (k = i; k <= l; k++) {
  1423.             putchar(itobase36(category[k]));
  1424.             if (((k % 10) == 0) && (k < l)) putchar(' ');
  1425.             }
  1426.           putchar('\n');
  1427.           }
  1428.  
  1429.         for (j = 1; j <= tr->mxtips; j++) {
  1430.           printf("%s   ", tr->nodep[j]->name);
  1431.           for (k = i; k <= l; k++) {
  1432.             ch = y[j][k];
  1433.             if ((j > 1) && (ch == y[1][k])) ch = '.';
  1434.             putchar(ch);
  1435.             if (((k % 10) == 0) && (k < l)) putchar(' ');
  1436.             }
  1437.           putchar('\n');
  1438.           }
  1439.         putchar('\n');
  1440.         }
  1441.       }
  1442.  
  1443.     for (j = 1; j <= tr->mxtips; j++)    /* Convert characters to meanings */
  1444.       for (i = 1; i <= sites; i++)  y[j][i] = meaning[y[j][i]];
  1445.  
  1446.   } /* getdata */
  1447.  
  1448.  
  1449. void getinput (tr)
  1450.     tree  *tr;
  1451.   { /* getinput */
  1452.     getnums(tr);                      if (anerror) return;
  1453.     getoptions(tr);                   if (anerror) return;
  1454.     if (! freqsfrom) getbasefreqs();  if (anerror) return;
  1455.     getyspace();                      if (anerror) return;
  1456.     setupTree(tr);                    if (anerror) return;
  1457.     getdata(tr);                      if (anerror) return;
  1458.   } /* getinput */
  1459.  
  1460.  
  1461. void sitesort ()
  1462.     /* Shell sort keeping sites, weights in same order */
  1463.   { /* sitesort */
  1464.     int  gap, i, j, jj, jg, k;
  1465.     boolean  flip, tied;
  1466.  
  1467.     for (gap = sites/2; gap > 0; gap /= 2) {
  1468.       for (i = gap + 1; i <= sites; i++) {
  1469.         j = i - gap;
  1470.  
  1471.         do {
  1472.           jj = alias[j];
  1473.           jg = alias[j+gap];
  1474.           flip = (category[jj] >  category[jg]);
  1475.           tied = (category[jj] == category[jg]);
  1476.           for (k = 1; (k <= numsp) && tied; k++) {
  1477.             flip = (y[k][jj] >  y[k][jg]);
  1478.             tied = (y[k][jj] == y[k][jg]);
  1479.             }
  1480.           if (flip) {
  1481.             alias[j]     = jg;
  1482.             alias[j+gap] = jj;
  1483.             j -= gap;
  1484.             }
  1485.           } while (flip && (j > 0));
  1486.  
  1487.         }  /* for (i ...   */
  1488.       }    /* for (gap ... */
  1489.   } /* sitesort */
  1490.  
  1491.  
  1492. void sitecombcrunch ()
  1493.     /* combine sites that have identical patterns (and nonzero weight) */
  1494.   { /* sitecombcrunch */
  1495.     int  i, sitei, j, sitej, k;
  1496.     boolean  tied;
  1497.  
  1498.     i = 0;
  1499.     alias[0] = alias[1];
  1500.     aliasweight[0] = 0;
  1501.  
  1502.     for (j = 1; j <= sites; j++) {
  1503.       sitei = alias[i];
  1504.       sitej = alias[j];
  1505.       tied = (category[sitei] == category[sitej]);
  1506.  
  1507.       for (k = 1; tied && (k <= numsp); k++)
  1508.         tied = (y[k][sitei] == y[k][sitej]);
  1509.  
  1510.       if (tied) {
  1511.         aliasweight[i] += weight[sitej];
  1512.         }
  1513.       else {
  1514.         if (aliasweight[i] > 0) i++;
  1515.         aliasweight[i] = weight[sitej];
  1516.         alias[i] = sitej;
  1517.         }
  1518.       }
  1519.  
  1520.     endsite = i;
  1521.     if (aliasweight[i] > 0) endsite++;
  1522.   } /* sitecombcrunch */
  1523.  
  1524.  
  1525. void makeweights ()
  1526.     /* make up weights vector to avoid duplicate computations */
  1527.   { /* makeweights */
  1528.     int  i;
  1529.  
  1530.     for (i = 1; i <= sites; i++)  alias[i] = i;
  1531.     sitesort();
  1532.     sitecombcrunch();
  1533.     if (endsite > maxpatterns) {
  1534.       printf("ERROR: Too many patterns in data\n");
  1535.       printf("       Increase maxpatterns to at least %d\n", endsite);
  1536.       anerror = TRUE;
  1537.       }
  1538.     else {
  1539.       printf("Analyzing %d distinct data patterns (columns)\n\n", endsite);
  1540.       }
  1541.   } /* makeweights */
  1542.  
  1543.  
  1544. void makevalues (tr)
  1545.     tree  *tr;
  1546.     /* set up fractional likelihoods at tips */
  1547.   { /* makevalues */
  1548.     double  temp;
  1549.     int  i, j, k;
  1550.  
  1551.     for (i = 1; i <= tr->mxtips; i++) {    /* Pack and move tip data */
  1552.       for (j = 0; j < endsite; j++)  y[i-1][j] = y[i][alias[j]];
  1553.       tr->nodep[i]->tip = &(y[i-1][0]);
  1554.       }
  1555.  
  1556.     totalwrate = 0.0;
  1557.     for (k = 0; k < endsite; k++) {
  1558.       catnumb[k] = i = category[alias[k]];
  1559.       ratvalue[k] = temp = rate[i];
  1560.       totalwrate += wgt_rate[k] = temp * aliasweight[k];
  1561.       wgt_rate2[k] = temp * temp * aliasweight[k];
  1562.       }
  1563.   } /* makevalues */
  1564.  
  1565.  
  1566. void empiricalfreqs (tr)
  1567.     tree  *tr;
  1568.     /* Get empirical base frequencies from the data */
  1569.   { /* empiricalfreqs */
  1570.     double  sum, suma, sumc, sumg, sumt, wj, fa, fc, fg, ft;
  1571.     int  i, j, k, code;
  1572.     char *yptr;
  1573.  
  1574.     freqa = 0.25;
  1575.     freqc = 0.25;
  1576.     freqg = 0.25;
  1577.     freqt = 0.25;
  1578.     for (k = 1; k <= 8; k++) {
  1579.       suma = 0.0;
  1580.       sumc = 0.0;
  1581.       sumg = 0.0;
  1582.       sumt = 0.0;
  1583.       for (i = 1; i <= tr->mxtips; i++) {
  1584.         yptr = tr->nodep[i]->tip;
  1585.         for (j = 0; j < endsite; j++) {
  1586.           code = *yptr++;
  1587.           fa = freqa * ( code       & 1);
  1588.           fc = freqc * ((code >> 1) & 1);
  1589.           fg = freqg * ((code >> 2) & 1);
  1590.           ft = freqt * ((code >> 3) & 1);
  1591.           wj = aliasweight[j] / (fa + fc + fg + ft);
  1592.           suma += wj * fa;
  1593.           sumc += wj * fc;
  1594.           sumg += wj * fg;
  1595.           sumt += wj * ft;
  1596.           }
  1597.         }
  1598.       sum = suma + sumc + sumg + sumt;
  1599.       freqa = suma / sum;
  1600.       freqc = sumc / sum;
  1601.       freqg = sumg / sum;
  1602.       freqt = sumt / sum;
  1603.       }
  1604.   } /* empiricalfreqs */
  1605.  
  1606.  
  1607. xarray *setupxarray ()
  1608.   { /* setupxarray */
  1609.     xarray  *x;
  1610.     xtype  *data;
  1611.  
  1612.     x = (xarray *) malloc((unsigned long) sizeof(xarray));
  1613.     if (x) {
  1614.       data = (xtype *) malloc((unsigned long) (4 * endsite * sizeof(xtype)));
  1615.       if (data) {
  1616.         x->a = data;
  1617.         x->c = data += endsite;
  1618.         x->g = data += endsite;
  1619.         x->t = data +  endsite;
  1620.         x->prev = x->next = x;
  1621.         x->owner = (node *) NULL;
  1622.         }
  1623.       else {
  1624.         free((char *) x);
  1625.         return (xarray *) NULL;
  1626.         }
  1627.       }
  1628.     return x;
  1629.   } /* setupxarray */
  1630.  
  1631.  
  1632. void linkxarray (req, min, freexptr, usedxptr)
  1633.     int  req, min;
  1634.     xarray **freexptr, **usedxptr;
  1635.     /*  Link a set of xarrays */
  1636.   { /* linkxarray */
  1637.     xarray  *first, *prev, *x;
  1638.     int  i;
  1639.  
  1640.     first = prev = (xarray *) NULL;
  1641.     i = 0;
  1642.  
  1643.     do {
  1644.       x = setupxarray();
  1645.       if (x) {
  1646.         if (! first) first = x;
  1647.         else {
  1648.           prev->next = x;
  1649.           x->prev = prev;
  1650.           }
  1651.         prev = x;
  1652.         i++;
  1653.         }
  1654.       else {
  1655.         printf("ERROR: Failure to get requested xarray memory\n");
  1656.         if (i < min) {
  1657.           anerror = TRUE;
  1658.           return;
  1659.           }
  1660.         }
  1661.       } while ((i < req) && x);
  1662.  
  1663.     if (first) {
  1664.       first->prev = prev;
  1665.       prev->next = first;
  1666.       }
  1667.  
  1668.     *freexptr = first;
  1669.     *usedxptr = (xarray *) NULL;
  1670.   } /* linkxarray */
  1671.  
  1672.  
  1673. void setupnodex (tr)
  1674.     tree  *tr;
  1675.   { /* setupnodex */
  1676.     nodeptr  p;
  1677.     int  i;
  1678.  
  1679.     for (i = tr->mxtips + 1; (i <= 2*(tr->mxtips) - 2); i++) {
  1680.       p = tr->nodep[i];
  1681.       if (anerror = !(p->x = setupxarray())) {
  1682.         printf("ERROR: Failure to get internal node xarray memory\n");
  1683.         return;
  1684.         }
  1685.       }
  1686.   } /* setupnodex */
  1687.  
  1688.  
  1689. xarray *getxtip (p)
  1690.     nodeptr  p;
  1691.   { /* getxtip */
  1692.     xarray  *new;
  1693.     boolean  splice;
  1694.  
  1695.     if (! p) return (xarray *) NULL;
  1696.  
  1697.     splice = FALSE;
  1698.  
  1699.     if (p->x) {                  /* array is there; move to tail of list */
  1700.       new = p->x;
  1701.       if (new == new->prev) ;             /* linked to self; leave it */
  1702.       else if (new == usedxtip) usedxtip = usedxtip->next; /* at head */
  1703.       else if (new == usedxtip->prev) ;   /* already at tail */
  1704.       else {                              /* move to tail of list */
  1705.         new->prev->next = new->next;
  1706.         new->next->prev = new->prev;
  1707.         splice = TRUE;
  1708.         }
  1709.       }
  1710.  
  1711.     else if (freextip) {                 /* take from unused list */
  1712.       p->x = new = freextip;
  1713.       new->owner = p;
  1714.       if (new->prev != new) {            /* not only member of freelist */
  1715.         new->prev->next = new->next;
  1716.         new->next->prev = new->prev;
  1717.         freextip = new->next;
  1718.         }
  1719.       else
  1720.         freextip = (xarray *) NULL;
  1721.  
  1722.       splice = TRUE;
  1723.       }
  1724.  
  1725.     else if (usedxtip) {                 /* take from head of used list */
  1726.       usedxtip->owner->x = (xarray *) NULL;
  1727.       p->x = new = usedxtip;
  1728.       new->owner = p;
  1729.       usedxtip = usedxtip->next;
  1730.       }
  1731.  
  1732.     else {
  1733.       printf("ERROR: Unable to locate memory for tip %d.\n", p->number);
  1734.       anerror = TRUE;
  1735.       exit(1);
  1736.       }
  1737.  
  1738.     if (splice) {
  1739.       if (usedxtip) {                  /* list is not empty */
  1740.         usedxtip->prev->next = new;
  1741.         new->prev = usedxtip->prev;
  1742.         usedxtip->prev = new;
  1743.         new->next = usedxtip;
  1744.         }
  1745.       else
  1746.         usedxtip = new->prev = new->next = new;
  1747.       }
  1748.  
  1749.     return  new;
  1750.   } /* getxtip */
  1751.  
  1752.  
  1753. xarray *getxnode (p)
  1754.     nodeptr  p;
  1755.     /* Ensure that internal node p has memory */
  1756.   { /* getxnode */
  1757.     nodeptr  s;
  1758.  
  1759.     if (! (p->x)) {  /*  Move likelihood array on this node to sector p */
  1760.       if ((s = p->next)->x || (s = s->next)->x) {
  1761.         p->x = s->x;
  1762.         s->x = (xarray *) NULL;
  1763.         }
  1764.       else {
  1765.         printf("ERROR: Unable to locate memory at node %d.\n", p->number);
  1766.         exit(1);
  1767.         }
  1768.       }
  1769.     return  p->x;
  1770.   } /* getxnode */
  1771.  
  1772.  
  1773. void newview (p)                      /*  Update likelihoods at node */
  1774.     nodeptr  p;
  1775.   { /* newview */
  1776.     double   z1, lz1, xvlz1, z2, lz2, xvlz2;
  1777.     nodeptr  q, r;
  1778.     xtype   *x1a, *x1c, *x1g, *x1t, *x2a, *x2c, *x2g, *x2t,
  1779.             *x3a, *x3c, *x3g, *x3t;
  1780.     int  i;
  1781.  
  1782.     if (p->tip) {             /*  Make sure that data are at tip */
  1783.       int code;
  1784.       char *yptr;
  1785.  
  1786.       if (p->x) return;       /*  They are already there */
  1787.  
  1788.       (void) getxtip(p);      /*  They are not, so get memory */
  1789.       x3a = &(p->x->a[0]);    /*  Move tip data to xarray */
  1790.       x3c = &(p->x->c[0]);
  1791.       x3g = &(p->x->g[0]);
  1792.       x3t = &(p->x->t[0]);
  1793.       yptr = p->tip;
  1794.       for (i = 0; i < endsite; i++) {
  1795.         code = *yptr++;
  1796.         *x3a++ =  code       & 1;
  1797.         *x3c++ = (code >> 1) & 1;
  1798.         *x3g++ = (code >> 2) & 1;
  1799.         *x3t++ = (code >> 3) & 1;
  1800.         }
  1801.       return;
  1802.       }
  1803.  
  1804. /*  Internal node needs update */
  1805.  
  1806.     q = p->next->back;
  1807.     r = p->next->next->back;
  1808.  
  1809.     while ((! p->x) || (! q->x) || (! r->x)) {
  1810.       if (! q->x) newview(q);
  1811.       if (! r->x) newview(r);
  1812.       if (! p->x) (void) getxnode(p);
  1813.       }
  1814.  
  1815.     x1a = &(q->x->a[0]);
  1816.     x1c = &(q->x->c[0]);
  1817.     x1g = &(q->x->g[0]);
  1818.     x1t = &(q->x->t[0]);
  1819.     z1 = q->z;
  1820.     lz1 = (z1 > zmin) ? log(z1) : log(zmin);
  1821.     xvlz1 = xv * lz1;
  1822.  
  1823.     x2a = &(r->x->a[0]);
  1824.     x2c = &(r->x->c[0]);
  1825.     x2g = &(r->x->g[0]);
  1826.     x2t = &(r->x->t[0]);
  1827.     z2 = r->z;
  1828.     lz2 = (z2 > zmin) ? log(z2) : log(zmin);
  1829.     xvlz2 = xv * lz2;
  1830.  
  1831.     x3a = &(p->x->a[0]);
  1832.     x3c = &(p->x->c[0]);
  1833.     x3g = &(p->x->g[0]);
  1834.     x3t = &(p->x->t[0]);
  1835.  
  1836.     { double  zz1table[maxcategories+1], zv1table[maxcategories+1],
  1837.               zz2table[maxcategories+1], zv2table[maxcategories+1],
  1838.               *zz1ptr, *zv1ptr, *zz2ptr, *zv2ptr, *rptr;
  1839.       double  fx1r, fx1y, fx1n, suma1, sumg1, sumc1, sumt1,
  1840.               fx2r, fx2y, fx2n, ki, tempi, tempj;
  1841.       int  *cptr;
  1842.  
  1843. #     ifdef Vectorize
  1844.         double zz1[maxsites], zv1[maxsites], zz2[maxsites], zv2[maxsites];
  1845.         int  cat;
  1846. #     else
  1847.         double zz1, zv1, zz2, zv2;
  1848.         int cat;
  1849. #     endif
  1850.  
  1851.       rptr   = &(rate[1]);
  1852.       zz1ptr = &(zz1table[1]);
  1853.       zv1ptr = &(zv1table[1]);
  1854.       zz2ptr = &(zz2table[1]);
  1855.       zv2ptr = &(zv2table[1]);
  1856.  
  1857. #     ifdef Vectorize
  1858. #       pragma IVDEP
  1859. #     endif
  1860.  
  1861.       for (i = 1; i <= categs; i++) {   /* exps for each category */
  1862.         ki = *rptr++;
  1863.         *zz1ptr++ = exp(ki *   lz1);
  1864.         *zv1ptr++ = exp(ki * xvlz1);
  1865.         *zz2ptr++ = exp(ki *   lz2);
  1866.         *zv2ptr++ = exp(ki * xvlz2);
  1867.         }
  1868.  
  1869.       cptr = &(catnumb[0]);
  1870.  
  1871. #     ifdef Vectorize
  1872. #       pragma IVDEP
  1873.         for (i = 0; i < endsite; i++) {
  1874.           cat = *cptr++;
  1875.           zz1[i] = zz1table[cat];
  1876.           zv1[i] = zv1table[cat];
  1877.           zz2[i] = zz2table[cat];
  1878.           zv2[i] = zv2table[cat];
  1879.         }
  1880.  
  1881. #       pragma IVDEP
  1882.         for (i = 0; i < endsite; i++) {
  1883.           fx1r = freqa * *x1a + freqg * *x1g;
  1884.           fx1y = freqc * *x1c + freqt * *x1t;
  1885.           fx1n = fx1r + fx1y;
  1886.           tempi = fx1r * invfreqr;
  1887.           tempj = zv1[i] * (tempi-fx1n) + fx1n;
  1888.           suma1 = zz1[i] * (*x1a++ - tempi) + tempj;
  1889.           sumg1 = zz1[i] * (*x1g++ - tempi) + tempj;
  1890.           tempi = fx1y * invfreqy;
  1891.           tempj = zv1[i] * (tempi-fx1n) + fx1n;
  1892.           sumc1 = zz1[i] * (*x1c++ - tempi) + tempj;
  1893.           sumt1 = zz1[i] * (*x1t++ - tempi) + tempj;
  1894.  
  1895.           fx2r = freqa * *x2a + freqg * *x2g;
  1896.           fx2y = freqc * *x2c + freqt * *x2t;
  1897.           fx2n = fx2r + fx2y;
  1898.           tempi = fx2r * invfreqr;
  1899.           tempj = zv2[i] * (tempi-fx2n) + fx2n;
  1900.           *x3a++ = suma1 * (zz2[i] * (*x2a++ - tempi) + tempj);
  1901.           *x3g++ = sumg1 * (zz2[i] * (*x2g++ - tempi) + tempj);
  1902.           tempi = fx2y * invfreqy;
  1903.           tempj = zv2[i] * (tempi-fx2n) + fx2n;
  1904.           *x3c++ = sumc1 * (zz2[i] * (*x2c++ - tempi) + tempj);
  1905.           *x3t++ = sumt1 * (zz2[i] * (*x2t++ - tempi) + tempj);
  1906.           }
  1907.  
  1908. #     else  /*  Not Vectorize  */
  1909.         for (i = 0; i < endsite; i++) {
  1910.           cat = *cptr++;
  1911.           zz1 = zz1table[cat];
  1912.           zv1 = zv1table[cat];
  1913.           fx1r = freqa * *x1a + freqg * *x1g;
  1914.           fx1y = freqc * *x1c + freqt * *x1t;
  1915.           fx1n = fx1r + fx1y;
  1916.           tempi = fx1r * invfreqr;
  1917.           tempj = zv1 * (tempi-fx1n) + fx1n;
  1918.           suma1 = zz1 * (*x1a++ - tempi) + tempj;
  1919.           sumg1 = zz1 * (*x1g++ - tempi) + tempj;
  1920.           tempi = fx1y * invfreqy;
  1921.           tempj = zv1 * (tempi-fx1n) + fx1n;
  1922.           sumc1 = zz1 * (*x1c++ - tempi) + tempj;
  1923.           sumt1 = zz1 * (*x1t++ - tempi) + tempj;
  1924.  
  1925.           zz2 = zz2table[cat];
  1926.           zv2 = zv2table[cat];
  1927.           fx2r = freqa * *x2a + freqg * *x2g;
  1928.           fx2y = freqc * *x2c + freqt * *x2t;
  1929.           fx2n = fx2r + fx2y;
  1930.           tempi = fx2r * invfreqr;
  1931.           tempj = zv2 * (tempi-fx2n) + fx2n;
  1932.           *x3a++ = suma1 * (zz2 * (*x2a++ - tempi) + tempj);
  1933.           *x3g++ = sumg1 * (zz2 * (*x2g++ - tempi) + tempj);
  1934.           tempi = fx2y * invfreqy;
  1935.           tempj = zv2 * (tempi-fx2n) + fx2n;
  1936.           *x3c++ = sumc1 * (zz2 * (*x2c++ - tempi) + tempj);
  1937.           *x3t++ = sumt1 * (zz2 * (*x2t++ - tempi) + tempj);
  1938.           }
  1939. #     endif  /*  Vectorize or not  */
  1940.  
  1941.       }
  1942.   } /* newview */
  1943.  
  1944.  
  1945. double evaluate (tr, p)
  1946.     tree  *tr;
  1947.     nodeptr  p;
  1948.   { /* evaluate */
  1949.     double   sum, z, lz, xvlz,
  1950.              ki, fx1a, fx1c, fx1g, fx1t, fx1r, fx1y, fx2r, fx2y,
  1951.              suma, sumb, sumc, term;
  1952.  
  1953. #   ifdef Vectorize
  1954.        double   zz[maxsites],zv[maxsites];
  1955. #   else
  1956.        double   zz,zv;
  1957. #   endif
  1958.  
  1959.     double   zztable[maxcategories+1], zvtable[maxcategories+1],
  1960.             *zzptr, *zvptr;
  1961.     double  *log_f, *rptr;
  1962.     xtype   *x1a, *x1c, *x1g, *x1t, *x2a, *x2c, *x2g, *x2t;
  1963.     nodeptr  q;
  1964.     int  cat, *cptr, i, *wptr;
  1965.  
  1966.     q = p->back;
  1967.     while ((! p->x) || (! q->x)) {
  1968.       if (! (p->x)) newview(p);
  1969.       if (! (q->x)) newview(q);
  1970.       }
  1971.  
  1972.     x1a = &(p->x->a[0]);
  1973.     x1c = &(p->x->c[0]);
  1974.     x1g = &(p->x->g[0]);
  1975.     x1t = &(p->x->t[0]);
  1976.  
  1977.     x2a = &(q->x->a[0]);
  1978.     x2c = &(q->x->c[0]);
  1979.     x2g = &(q->x->g[0]);
  1980.     x2t = &(q->x->t[0]);
  1981.  
  1982.     z = p->z;
  1983.     if (z < zmin) z = zmin;
  1984.     lz = log(z);
  1985.     xvlz = xv * lz;
  1986.  
  1987.     rptr  = &(rate[1]);
  1988.     zzptr = &(zztable[1]);
  1989.     zvptr = &(zvtable[1]);
  1990.  
  1991. #   ifdef Vectorize
  1992. #     pragma IVDEP
  1993. #   endif
  1994.  
  1995.     for (i = 1; i <= categs; i++) {
  1996.       ki = *rptr++;
  1997.       *zzptr++ = exp(ki *   lz);
  1998.       *zvptr++ = exp(ki * xvlz);
  1999.       }
  2000.  
  2001.     wptr = &(aliasweight[0]);
  2002.     cptr = &(catnumb[0]);
  2003.     log_f = tr->log_f;
  2004.     tr->log_f_valid = TRUE;
  2005.     sum = 0.0;
  2006.  
  2007. #   ifdef Vectorize
  2008. #     pragma IVDEP
  2009.       for (i = 0; i < endsite; i++) {
  2010.         cat   = *cptr++;
  2011.         zz[i] = zztable[cat];
  2012.         zv[i] = zvtable[cat];
  2013.       }
  2014.  
  2015. #     pragma IVDEP
  2016.       for (i = 0; i < endsite; i++) {
  2017.         fx1a = freqa * *x1a++;
  2018.         fx1g = freqg * *x1g++;
  2019.         fx1c = freqc * *x1c++;
  2020.         fx1t = freqt * *x1t++;
  2021.         suma = fx1a * *x2a + fx1c * *x2c + fx1g * *x2g + fx1t * *x2t;
  2022.         fx2r = freqa * *x2a++ + freqg * *x2g++;
  2023.         fx2y = freqc * *x2c++ + freqt * *x2t++;
  2024.         fx1r = fx1a + fx1g;
  2025.         fx1y = fx1c + fx1t;
  2026.         sumc = (fx1r + fx1y) * (fx2r + fx2y);
  2027.         sumb = fx1r * fx2r * invfreqr + fx1y * fx2y * invfreqy;
  2028.         suma -= sumb;
  2029.         sumb -= sumc;
  2030.         term = log(zz[i] * suma + zv[i] * sumb + sumc);
  2031.         sum += *wptr++ * term;
  2032.         *log_f++ = term;
  2033.         }
  2034.  
  2035. #   else  /*  Not Vectorize  */
  2036.       for (i = 0; i < endsite; i++) {
  2037.         cat  = *cptr++;
  2038.         zz   = zztable[cat];
  2039.         zv   = zvtable[cat];
  2040.         fx1a = freqa * *x1a++;
  2041.         fx1g = freqg * *x1g++;
  2042.         fx1c = freqc * *x1c++;
  2043.         fx1t = freqt * *x1t++;
  2044.         suma = fx1a * *x2a + fx1c * *x2c + fx1g * *x2g + fx1t * *x2t;
  2045.         fx2r = freqa * *x2a++ + freqg * *x2g++;
  2046.         fx2y = freqc * *x2c++ + freqt * *x2t++;
  2047.         fx1r = fx1a + fx1g;
  2048.         fx1y = fx1c + fx1t;
  2049.         sumc = (fx1r + fx1y) * (fx2r + fx2y);
  2050.         sumb = fx1r * fx2r * invfreqr + fx1y * fx2y * invfreqy;
  2051.         suma -= sumb;
  2052.         sumb -= sumc;
  2053.         term = log(zz * suma + zv * sumb + sumc);
  2054.         sum += *wptr++ * term;
  2055.         *log_f++ = term;
  2056.         }
  2057. #   endif  /*  Vectorize or not  */
  2058.  
  2059.     tr->likelihood = sum;
  2060.     return  sum;
  2061.   } /* evaluate */
  2062.  
  2063.  
  2064. #if 0                           /*  This is not currently used */
  2065. double evaluateslope (p, q, z)  /*  d(L)/d(lz) */
  2066.     nodeptr  p, q;
  2067.     double  z;
  2068.   { /* evaluateslope */
  2069.     double   dLdlz, lz, xvlz;
  2070.     double   sumai[maxpatterns], sumbi[maxpatterns], sumci[maxpatterns],
  2071.             *sumaptr, *sumbptr, *sumcptr;
  2072.     double   zztable[maxcategories+1], zvtable[maxcategories+1],
  2073.             *zzptr, *zvptr;
  2074.     double   ki, fx1a, fx1c, fx1g, fx1t, fx1r, fx1y, fx2r, fx2y,
  2075.              suma, sumb, sumc;
  2076.     double  *rptr, *wrptr;
  2077.     xtype   *x1a, *x1c, *x1g, *x1t, *x2a, *x2c, *x2g, *x2t;
  2078.     int  cat, *cptr, i;
  2079.  
  2080.  
  2081.     while ((! p->x) || (! q->x)) {
  2082.       if (! p->x) newview(p);
  2083.       if (! q->x) newview(q);
  2084.       }
  2085.  
  2086.     x1a = &(p->x->a[0]);
  2087.     x1c = &(p->x->c[0]);
  2088.     x1g = &(p->x->g[0]);
  2089.     x1t = &(p->x->t[0]);
  2090.  
  2091.     x2a = &(q->x->a[0]);
  2092.     x2c = &(q->x->c[0]);
  2093.     x2g = &(q->x->g[0]);
  2094.     x2t = &(q->x->t[0]);
  2095.  
  2096.     sumaptr = &(sumai[0]);
  2097.     sumbptr = &(sumbi[0]);
  2098.     sumcptr = &(sumci[0]);
  2099.  
  2100.     for (i = 0; i < endsite; i++) {
  2101.       fx1a = freqa * *x1a++;
  2102.       fx1g = freqg * *x1g++;
  2103.       fx1c = freqc * *x1c++;
  2104.       fx1t = freqt * *x1t++;
  2105.       suma   = fx1a * *x2a + fx1c * *x2c + fx1g * *x2g + fx1t * *x2t;
  2106.       fx2r = freqa * *x2a++ + freqg * *x2g++;
  2107.       fx2y = freqc * *x2c++ + freqt * *x2t++;
  2108.       fx1r = fx1a + fx1g;
  2109.       fx1y = fx1c + fx1t;
  2110.       *sumcptr++ = sumc = (fx1r + fx1y) * (fx2r + fx2y);
  2111.       sumb       = fx1r * fx2r * invfreqr + fx1y * fx2y * invfreqy;
  2112.       *sumaptr++ = suma - sumb;
  2113.       *sumbptr++ = sumb - sumc;
  2114.       }
  2115.  
  2116.     if (z < zmin) z = zmin;
  2117.     lz = log(z);
  2118.     xvlz = xv * lz;
  2119.  
  2120.     rptr  = &(rate[1]);
  2121.     zzptr = &(zztable[1]);
  2122.     zvptr = &(zvtable[1]);
  2123.     for (i = 1; i <= categs; i++) {
  2124.       ki = *rptr++;
  2125.       *zzptr++ = exp(ki *   lz);
  2126.       *zvptr++ = exp(ki * xvlz);
  2127.       }
  2128.  
  2129.     sumaptr = &(sumai[0]);
  2130.     sumbptr = &(sumbi[0]);
  2131.     sumcptr = &(sumci[0]);
  2132.     cptr  = &(catnumb[0]);
  2133.     wrptr = &(wgt_rate[0]);
  2134.     dLdlz = 0.0;
  2135.  
  2136.     for (i = 0; i < endsite; i++) {
  2137.       cat    = *cptr++;
  2138.       suma   = *sumaptr++ * zztable[cat];
  2139.       sumb   = *sumbptr++ * zvtable[cat];
  2140.       dLdlz += *wrptr++ * (suma + sumb*xv) / (suma + sumb + *sumcptr++);
  2141.       }
  2142.  
  2143.     return  dLdlz;
  2144.   } /* evaluateslope */
  2145. #endif
  2146.  
  2147.  
  2148. #ifdef  SmallMemorySystem 
  2149. static double   abi[maxpatterns], bci[maxpatterns], sumci[maxpatterns];
  2150. static double   zztable[maxcategories+1], zvtable[maxcategories+1];
  2151. #endif
  2152.  
  2153. double makenewz (p, q, z0, maxiter)
  2154.     nodeptr  p, q;
  2155.     double  z0;
  2156.     int  maxiter;
  2157.   { /* makenewz */
  2158. #ifndef SmallMemorySystem 
  2159.     double   abi[maxpatterns], bci[maxpatterns], sumci[maxpatterns];
  2160.     double   zztable[maxcategories+1], zvtable[maxcategories+1];
  2161. #endif
  2162.     double *abptr, *bcptr, *sumcptr;
  2163.     double   dlnLidlz, dlnLdlz, d2lnLdlz2, z, zprev, zstep, lz, xvlz,
  2164.              ki, suma, sumb, sumc, ab, bc, inv_Li, t1, t2,
  2165.              fx1a, fx1c, fx1g, fx1t, fx1r, fx1y, fx2r, fx2y;
  2166.     double    *zzptr, *zvptr;
  2167.     double  *rptr, *wrptr, *wr2ptr;
  2168.     xtype   *x1a, *x1c, *x1g, *x1t, *x2a, *x2c, *x2g, *x2t;
  2169.     int    cat, *cptr, i, curvatOK;
  2170.  
  2171.  
  2172.     while ((! p->x) || (! q->x)) {
  2173.       if (! (p->x)) newview(p);
  2174.       if (! (q->x)) newview(q);
  2175.       }
  2176.  
  2177.     x1a = &(p->x->a[0]);
  2178.     x1c = &(p->x->c[0]);
  2179.     x1g = &(p->x->g[0]);
  2180.     x1t = &(p->x->t[0]);
  2181.     x2a = &(q->x->a[0]);
  2182.     x2c = &(q->x->c[0]);
  2183.     x2g = &(q->x->g[0]);
  2184.     x2t = &(q->x->t[0]);
  2185.  
  2186.     abptr = &(abi[0]);
  2187.     bcptr = &(bci[0]);
  2188.     sumcptr = &(sumci[0]);
  2189.  
  2190. #   ifdef Vectorize
  2191. #     pragma IVDEP
  2192. #   endif
  2193.  
  2194.     for (i = 0; i < endsite; i++) {
  2195.       fx1a = freqa * *x1a++;
  2196.       fx1g = freqg * *x1g++;
  2197.       fx1c = freqc * *x1c++;
  2198.       fx1t = freqt * *x1t++;
  2199.       suma = fx1a * *x2a + fx1c * *x2c + fx1g * *x2g + fx1t * *x2t;
  2200.       fx2r = freqa * *x2a++ + freqg * *x2g++;
  2201.       fx2y = freqc * *x2c++ + freqt * *x2t++;
  2202.       fx1r = fx1a + fx1g;
  2203.       fx1y = fx1c + fx1t;
  2204.       *sumcptr++ = sumc = (fx1r + fx1y) * (fx2r + fx2y);
  2205.       sumb       = fx1r * fx2r * invfreqr + fx1y * fx2y * invfreqy;
  2206.       *abptr++   = suma - sumb;
  2207.       *bcptr++   = sumb - sumc;
  2208.       }
  2209.  
  2210.     z = z0;
  2211.     do {
  2212.       zprev = z;
  2213.       zstep = (1.0 - zmax) * z + zmin;
  2214.       curvatOK = FALSE;
  2215.  
  2216.       do {
  2217.         if (z < zmin) z = zmin;
  2218.         else if (z > zmax) z = zmax;
  2219.  
  2220.         lz    = log(z);
  2221.         xvlz  = xv * lz;
  2222.         rptr  = &(rate[1]);
  2223.         zzptr = &(zztable[1]);
  2224.         zvptr = &(zvtable[1]);
  2225.  
  2226. #       ifdef Vectorize
  2227. #         pragma IVDEP
  2228. #       endif
  2229.  
  2230.         for (i = 1; i <= categs; i++) {
  2231.           ki = *rptr++;
  2232.           *zzptr++ = exp(ki *   lz);
  2233.           *zvptr++ = exp(ki * xvlz);
  2234.           }
  2235.  
  2236.         abptr   = &(abi[0]);
  2237.         bcptr   = &(bci[0]);
  2238.         sumcptr = &(sumci[0]);
  2239.         cptr    = &(catnumb[0]);
  2240.         wrptr   = &(wgt_rate[0]);
  2241.         wr2ptr  = &(wgt_rate2[0]);
  2242.         dlnLdlz = 0.0;                 /*  = d(ln(likelihood))/d(lz) */
  2243.         d2lnLdlz2 = 0.0;               /*  = d2(ln(likelihood))/d(lz)2 */
  2244.  
  2245. #       ifdef Vectorize
  2246. #         pragma IVDEP
  2247. #       endif
  2248.  
  2249.         for (i = 0; i < endsite; i++) {
  2250.           cat    = *cptr++;              /*  ratecategory(i) */
  2251.           ab     = *abptr++ * zztable[cat];
  2252.           bc     = *bcptr++ * zvtable[cat];
  2253.           sumc   = *sumcptr++;
  2254.           inv_Li = 1.0/(ab + bc + sumc);
  2255.           t1     = ab * inv_Li;
  2256.           t2     = xv * bc * inv_Li;
  2257.           dlnLidlz   = t1 + t2;
  2258.           dlnLdlz   += *wrptr++  * dlnLidlz;
  2259.           d2lnLdlz2 += *wr2ptr++ * (t1 + xv * t2 - dlnLidlz * dlnLidlz);
  2260.           }
  2261.  
  2262.         if ((d2lnLdlz2 >= 0.0) && (z < zmax))
  2263.           zprev = z = 0.37 * z + 0.63;  /*  Bad curvature, shorten branch */
  2264.         else
  2265.           curvatOK = TRUE;
  2266.  
  2267.         } while (! curvatOK);
  2268.  
  2269.       if (d2lnLdlz2 < 0.0) {
  2270.         z *= exp(-dlnLdlz / d2lnLdlz2);
  2271.         if (z < zmin) z = zmin;
  2272.         if (z > 0.25 * zprev + 0.75)    /*  Limit steps toward z = 1.0 */
  2273.           z = 0.25 * zprev + 0.75;
  2274.         }
  2275.  
  2276.       if (z > zmax) z = zmax;
  2277.  
  2278.       } while ((--maxiter > 0) && (ABS(z - zprev) > zstep));
  2279.  
  2280.     return  z;
  2281.   } /* makenewz */
  2282.  
  2283.  
  2284. void update (tr, p)
  2285.     tree    *tr;
  2286.     nodeptr  p;
  2287.   { /* update */
  2288.     nodeptr  q;
  2289.     double   z0, z;
  2290.  
  2291.     q = p->back;
  2292.     z0 = q->z;
  2293.     p->z = q->z = z = makenewz(p, q, z0, newzpercycle);
  2294.     if (ABS(z - z0) > deltaz)  tr->smoothed = FALSE;
  2295.   } /* update */
  2296.  
  2297.  
  2298. void smooth (tr, p)
  2299.     tree    *tr;
  2300.     nodeptr  p;
  2301.   { /* smooth */
  2302.     update(tr, p);                       /*  Adjust branch */
  2303.     if (! p->tip) {                      /*  Adjust "descendents" */
  2304.       smooth(tr, p->next->back);
  2305.       smooth(tr, p->next->next->back);
  2306.  
  2307. #     if ReturnSmoothedView
  2308.         newview(p);
  2309. #     endif
  2310.       }
  2311.   } /* smooth */
  2312.  
  2313.  
  2314. void smoothTree (tr, maxtimes)
  2315.     tree    *tr;
  2316.     int  maxtimes;
  2317.   { /* smoothTree */
  2318.     nodeptr  p;
  2319.  
  2320.     p = tr->start;
  2321.  
  2322.     while (--maxtimes >= 0 && ! anerror) {
  2323.       tr->smoothed = TRUE;
  2324.       smooth(tr, p->back);
  2325.       if (! p->tip) {
  2326.         smooth(tr, p->next->back);
  2327.         smooth(tr, p->next->next->back);
  2328.         }
  2329.       if (tr->smoothed)  break;
  2330.       }
  2331.   } /* smoothTree */
  2332.  
  2333.  
  2334. void localSmooth (tr, p, maxtimes)    /* Smooth branches around p */
  2335.     tree    *tr;
  2336.     nodeptr  p;
  2337.     int  maxtimes;
  2338.   { /* localSmooth */
  2339.     nodeptr  pn, pnn;
  2340.  
  2341.     if (p->tip) return;               /* Should actually be an error */
  2342.  
  2343.     pn  = p->next;
  2344.     pnn = pn->next;
  2345.     while (--maxtimes >= 0) {
  2346.       tr->smoothed = TRUE;
  2347.       update(tr, p);     if (anerror) break;
  2348.       update(tr, pn);    if (anerror) break;
  2349.       update(tr, pnn);   if (anerror) break;
  2350.       if (tr->smoothed)  break;
  2351.       }
  2352.     tr->smoothed = FALSE;             /* Only smooth locally */
  2353.   } /* localSmooth */
  2354.  
  2355.  
  2356. void hookup (p, q, z)
  2357.     nodeptr  p, q;
  2358.     double   z;
  2359.   { /* hookup */
  2360.     p->back = q;
  2361.     q->back = p;
  2362.     p->z = q->z = z;
  2363.   } /* hookup */
  2364.  
  2365.  
  2366. void insert (tr, p, q, glob)   /* Insert node p into branch q <-> q->back */
  2367.     tree    *tr;
  2368.     nodeptr  p, q;
  2369.     boolean  glob;             /* Smooth tree globally? */
  2370.  
  2371. /*                q
  2372.                  /.
  2373.              add/ .
  2374.                /  .
  2375.              pn   .
  2376.     s ---- p      .remove
  2377.              pnn  .
  2378.                \  .
  2379.              add\ .
  2380.                  \.      pn  = p->next;
  2381.                   r      pnn = p->next->next;
  2382.  */
  2383.  
  2384.   { /* insert */
  2385.     nodeptr  r, s;
  2386.  
  2387.     r = q->back;
  2388.     s = p->back;
  2389.  
  2390. #   if BestInsertAverage && ! Master
  2391.     { double  zqr, zqs, zrs, lzqr, lzqs, lzrs, lzsum, lzq, lzr, lzs, lzmax;
  2392.  
  2393.       zqr = makenewz(q, r, q->z,     iterations);
  2394.       zqs = makenewz(q, s, defaultz, iterations);
  2395.       zrs = makenewz(r, s, defaultz, iterations);
  2396.  
  2397.       lzqr = (zqr > zmin) ? log(zqr) : log(zmin);  /* long branches */
  2398.       lzqs = (zqs > zmin) ? log(zqs) : log(zmin);
  2399.       lzrs = (zrs > zmin) ? log(zrs) : log(zmin);
  2400.       lzsum = 0.5 * (lzqr + lzqs + lzrs);
  2401.       lzq = lzsum - lzrs;
  2402.       lzr = lzsum - lzqs;
  2403.       lzs = lzsum - lzqr;
  2404.       lzmax = log(zmax);
  2405.       if      (lzq > lzmax) {lzq = lzmax; lzr = lzqr; lzs = lzqs;} /* short */
  2406.       else if (lzr > lzmax) {lzr = lzmax; lzq = lzqr; lzs = lzrs;}
  2407.       else if (lzs > lzmax) {lzs = lzmax; lzq = lzqs; lzr = lzrs;}
  2408.  
  2409.       hookup(p->next,       q, exp(lzq));
  2410.       hookup(p->next->next, r, exp(lzr));
  2411.       hookup(p,             s, exp(lzs));
  2412.       }
  2413.  
  2414. #   else
  2415.     { double  z;
  2416.       z = sqrt(q->z);
  2417.       hookup(p->next,       q, z);
  2418.       hookup(p->next->next, r, z);
  2419.       }
  2420.  
  2421. #   endif
  2422.  
  2423.     newview(p);         /*  Required so that sector p is valid at update */
  2424.     tr->opt_level = 0;
  2425.  
  2426. #   if ! Master         /*  Smoothings are done by slave */
  2427.       if (! glob)  localSmooth(tr, p, smoothings);  /* Smooth locale of p */
  2428.       if (glob)    smoothTree(tr, smoothings);      /* Smooth whole tree */
  2429.  
  2430. #   else
  2431.       tr->likelihood = unlikely;
  2432. #   endif
  2433.  
  2434.   } /* insert */
  2435.  
  2436.  
  2437. nodeptr  removeNode (tr, p)
  2438.     tree    *tr;
  2439.     nodeptr  p;
  2440.  
  2441. /*                q
  2442.                  .|
  2443.           remove. |
  2444.                .  |
  2445.              pn   |
  2446.     s ---- p      |add
  2447.              pnn  |
  2448.                .  |
  2449.           remove. |
  2450.                  .|      pn  = p->next;
  2451.                   r      pnn = p->next->next;
  2452.  */
  2453.  
  2454.     /* remove p and return where it was */
  2455.   { /* removeNode */
  2456.     double   zqr;
  2457.     nodeptr  q, r;
  2458.  
  2459.     q = p->next->back;
  2460.     r = p->next->next->back;
  2461.     zqr = q->z * r->z;
  2462. #   if ! Master
  2463.       hookup(q, r, makenewz(q, r, zqr, iterations));
  2464. #   else
  2465.       hookup(q, r, zqr);
  2466. #   endif
  2467.  
  2468.     p->next->next->back = p->next->back = (node *) NULL;
  2469.     return  q;
  2470.   } /* removeNode */
  2471.  
  2472.  
  2473. void initrav (p)
  2474.     nodeptr  p;
  2475.   { /* initrav */
  2476.     if (! p->tip) {
  2477.       initrav(p->next->back);
  2478.       initrav(p->next->next->back);
  2479.       newview(p);
  2480.       }
  2481.   } /* initrav */
  2482.  
  2483.  
  2484. nodeptr buildNewTip (tr, p)
  2485.     tree  *tr;
  2486.     nodeptr  p;
  2487.   { /* buildNewTip */
  2488.     nodeptr  q;
  2489.  
  2490.     q = tr->nodep[(tr->nextnode)++];
  2491.     hookup(p, q, defaultz);
  2492.     return  q;
  2493.   } /* buildNewTip */
  2494.  
  2495.  
  2496. void buildSimpleTree (tr, ip, iq, ir)
  2497.     tree  *tr;
  2498.     int    ip, iq, ir;
  2499.   { /* buildSimpleTree */
  2500.     /* p, q and r are tips meeting at s */
  2501.     nodeptr  p, s;
  2502.     int  i;
  2503.  
  2504.     i = MIN(ip, iq);
  2505.     if (ir < i)  i = ir; 
  2506.     tr->start = tr->nodep[i];
  2507.     tr->ntips = 3;
  2508.     p = tr->nodep[ip];
  2509.     hookup(p, tr->nodep[iq], defaultz);
  2510.     s = buildNewTip(tr, tr->nodep[ir]);
  2511.     insert(tr, s, p, FALSE);            /* Smoothing is local to s */
  2512.   } /* buildSimpleTree */
  2513.  
  2514.  
  2515. #ifndef HasStrLib
  2516.  
  2517. char * strchr (str, chr)
  2518.     char *str;
  2519.     int   chr;
  2520.  { /* strchr */
  2521.     int  c;
  2522.  
  2523.     while (c = *str)  {if (c == chr) return str; str++;}
  2524.     return  (char *) NULL;
  2525.  } /* strchr */
  2526.  
  2527.  
  2528. char * strstr (str1, str2)
  2529.     char *str1, *str2;
  2530.  { /* strstr */
  2531.     char *s1, *s2;
  2532.     int  c;
  2533.  
  2534.     while (*(s1 = str1)) {
  2535.       s2 = str2;
  2536.       do {
  2537.         if (! (c = *s2++))  return str1;
  2538.         } 
  2539.         while (*s1++ == c);
  2540.       str1++;
  2541.       }
  2542.     return  (char *) NULL;
  2543.  } /* strstr */
  2544.  
  2545. #endif
  2546.  
  2547. boolean readKeyValue (string, key, format, value)
  2548.     char *string, *key, *format;
  2549.     void *value;
  2550.   { /* readKeyValue */
  2551.  
  2552.     if (! (string = strstr(string, key)))  return FALSE;
  2553.     string += strlen(key);
  2554.     if (! (string = strchr(string, '=')))  return FALSE;
  2555.     string++;
  2556.     return  sscanf(string, format, value);  /* 1 if read, otherwise 0 */
  2557.   } /* readKeyValue */
  2558.  
  2559.  
  2560. #if Master || Slave
  2561.  
  2562. double str_readTreeLikelihood (treestr)
  2563.     char *treestr;
  2564.   { /* str_readTreeLikelihood */
  2565.     double lk1;
  2566.     char    *com, *com_end;
  2567.     boolean  readKeyValue();
  2568.  
  2569.     if ((com = strchr(treestr, '[')) && (com < strchr(treestr, '('))
  2570.                                      && (com_end = strchr(com, ']'))) {
  2571.       com++;
  2572.       *com_end = 0;
  2573.       if (readKeyValue(com, likelihood_key, "%lg", (void *) &(lk1))) {
  2574.         *com_end = ']';
  2575.         return lk1;
  2576.         }
  2577.       }
  2578.  
  2579.     fprintf(stderr, "ERROR reading likelihood in receiveTree\n");
  2580.     anerror = TRUE;
  2581.     return 1.0;
  2582.   } /* str_readTreeLikelihood */
  2583.  
  2584.  
  2585. void sendTree (comm, tr)
  2586.     comm_block *comm;
  2587.     tree *tr;
  2588.   { /* sendTree */
  2589.     char treestr[maxsp*(nmlngth+24)+100];
  2590.     char *treeString();
  2591. #   if Master
  2592.       void sendTreeNum();
  2593. #   endif
  2594.  
  2595.     comm->done_flag = tr->likelihood > 0.0;
  2596.     if (comm->done_flag)
  2597.       write_comm_msg(comm, NULL);
  2598.  
  2599.     else {
  2600. #     if Master
  2601.         if (send_ahead >= MAX_SEND_AHEAD) {
  2602.           double new_likelihood;
  2603.           int  n_to_get;
  2604.  
  2605.           n_to_get = (send_ahead+1)/2;
  2606.           sendTreeNum(n_to_get);
  2607.           send_ahead -= n_to_get;
  2608.           read_comm_msg(&comm_slave, treestr);
  2609.           new_likelihood = str_readTreeLikelihood(treestr);
  2610.           if (! best_tr_recv || (new_likelihood > best_lk_recv)) {
  2611.             if (best_tr_recv)  free(best_tr_recv);
  2612.             best_tr_recv = malloc((unsigned long) (strlen(treestr) + 1));
  2613.             strcpy(best_tr_recv, treestr);
  2614.             best_lk_recv = new_likelihood;
  2615.             }
  2616.           }
  2617.         send_ahead++;
  2618. #     endif           /*  End #if Master  */
  2619.  
  2620.       REPORT_SEND_TREE;
  2621.       (void) treeString(treestr, tr, tr->start->back, 1);
  2622.       write_comm_msg(comm, treestr);
  2623.       }
  2624.   } /* sendTree */
  2625.  
  2626.  
  2627. void  receiveTree (comm, tr)
  2628.     comm_block  *comm;
  2629.     tree        *tr;
  2630.   { /* receiveTree */
  2631.     char treestr[maxsp*(nmlngth+24)+100];
  2632.     void str_treeReadLen();
  2633.  
  2634.     read_comm_msg(comm, treestr);
  2635.     if (comm->done_flag)
  2636.       tr->likelihood = 1.0;
  2637.  
  2638.     else {
  2639. #     if Master
  2640.         if (best_tr_recv) {
  2641.           if (str_readTreeLikelihood(treestr) < best_lk_recv) {
  2642.             strcpy(treestr, best_tr_recv);  /* Overwrite new tree with best */
  2643.             }
  2644.           free(best_tr_recv);
  2645.           best_tr_recv = NULL;
  2646.           }
  2647. #     endif           /*  End #if Master  */
  2648.  
  2649.       str_treeReadLen(treestr, tr);
  2650.       }
  2651.   } /* receiveTree */
  2652.  
  2653.  
  2654. void requestForWork ()
  2655.   { /* requestForWork */
  2656.     p4_send(DNAML_REQUEST, DNAML_DISPATCHER_ID, NULL, 0);
  2657.   } /* requestForWork */
  2658. #endif                  /* End #if Master || Slave  */
  2659.  
  2660.  
  2661. #if Master
  2662. void sendTreeNum(n_to_get)
  2663.     int n_to_get;
  2664.   { /* sendTreeNum */
  2665.     char scr[512];
  2666.  
  2667.     sprintf(scr, "%d", n_to_get);
  2668.     p4_send(DNAML_NUM_TREE, DNAML_MERGER_ID, scr, strlen(scr)+1);
  2669.   } /* sendTreeNum */
  2670.  
  2671.  
  2672. void  getReturnedTrees (tr, bt, n_tree_sent)
  2673.     tree     *tr;
  2674.     bestlist *bt;
  2675.     int n_tree_sent; /* number of trees sent to slaves */
  2676.   { /* getReturnedTrees */
  2677.     void sendTreeNum(), receiveTree();
  2678.  
  2679.     sendTreeNum(send_ahead);
  2680.     send_ahead = 0;
  2681.  
  2682.     receiveTree(&comm_slave, tr);
  2683.     tr->smoothed = TRUE;
  2684.     (void) saveBestTree(bt, tr);
  2685.   } /* getReturnedTrees */
  2686. #endif
  2687.  
  2688.  
  2689. void cacheZ (tr)
  2690.     tree  *tr;
  2691.   { /* cacheZ */
  2692.     nodeptr  p;
  2693.     int  nodes;
  2694.  
  2695.     nodes = tr->mxtips  +  3 * (tr->mxtips - 2);
  2696.     p = tr->nodep[1];
  2697.     while (nodes-- > 0) {p->z0 = p->z; p++;}
  2698.   } /* cacheZ */
  2699.  
  2700.  
  2701. void restoreZ (tr)
  2702.     tree  *tr;
  2703.   { /* restoreZ */
  2704.     nodeptr  p;
  2705.     int  nodes;
  2706.  
  2707.     nodes = tr->mxtips  +  3 * (tr->mxtips - 2);
  2708.     p = tr->nodep[1];
  2709.     while (nodes-- > 0) {p->z = p->z0; p++;}
  2710.   } /* restoreZ */
  2711.  
  2712.  
  2713. testInsert (tr, p, q, bt, fast)
  2714.     tree     *tr;
  2715.     nodeptr   p, q;
  2716.     bestlist *bt;
  2717.     boolean   fast;
  2718.   { /* testInsert */
  2719.     double  qz;
  2720.     nodeptr  r;
  2721.  
  2722.     r = q->back;             /* Save original connection */
  2723.     qz = q->z;
  2724.     insert(tr, p, q, ! fast);
  2725.  
  2726. #   if ! Master
  2727.       (void) evaluate(tr, fast ? p->next->next : tr->start);
  2728.       (void) saveBestTree(bt, tr);
  2729. #   else  /* Master */
  2730.       tr->likelihood = unlikely;
  2731.       sendTree(&comm_slave, tr);
  2732. #   endif
  2733.  
  2734.     /* remove p from this branch */
  2735.  
  2736.     hookup(q, r, qz);
  2737.     p->next->next->back = p->next->back = (nodeptr) NULL;
  2738.     if (! fast) {            /* With fast add, other values are still OK */
  2739.       restoreZ(tr);          /*   Restore branch lengths */
  2740. #     if ! Master            /*   Regenerate x values */
  2741.         initrav(p->back);
  2742.         initrav(q);
  2743.         initrav(r);
  2744. #     endif
  2745.       }
  2746.   } /* testInsert */
  2747.  
  2748.  
  2749. int addTraverse (tr, p, q, mintrav, maxtrav, bt, fast)
  2750.     tree     *tr;
  2751.     nodeptr   p, q;
  2752.     int       mintrav, maxtrav;
  2753.     bestlist *bt;
  2754.     boolean   fast;
  2755.   { /* addTraverse */
  2756.     int  tested;
  2757.  
  2758.     tested = 0;
  2759.     if (--mintrav <= 0) {           /* Moved minimum distance? */
  2760.       testInsert(tr, p, q, bt, fast);
  2761.       tested++;
  2762.       }
  2763.  
  2764.     if ((! q->tip) && (--maxtrav > 0)) {    /* Continue traverse? */
  2765.       tested += addTraverse(tr, p, q->next->back,
  2766.                             mintrav, maxtrav, bt, fast);
  2767.       tested += addTraverse(tr, p, q->next->next->back,
  2768.                             mintrav, maxtrav, bt, fast);
  2769.       }
  2770.  
  2771.     return tested;
  2772.   } /* addTraverse */
  2773.  
  2774.  
  2775. int  rearrange (tr, p, mintrav, maxtrav, bt)
  2776.     tree     *tr;
  2777.     nodeptr   p;
  2778.     int       mintrav, maxtrav;
  2779.     bestlist *bt;
  2780.     /* rearranges the tree, globally or locally */
  2781.   { /* rearrange */
  2782.     double   p1z, p2z, q1z, q2z;
  2783.     nodeptr  p1, p2, q, q1, q2;
  2784.     int      tested, mintrav2;
  2785.  
  2786.     tested = 0;
  2787.     if (maxtrav < 1 || mintrav > maxtrav)  return tested;
  2788.  
  2789. /* Moving subtree forward in tree. */
  2790.  
  2791.     if (! p->tip) {
  2792.       p1 = p->next->back;
  2793.       p2 = p->next->next->back;
  2794.       if (! p1->tip || ! p2->tip) {
  2795.         p1z = p1->z;
  2796.         p2z = p2->z;
  2797.         (void) removeNode(tr, p);
  2798.         cacheZ(tr);
  2799.         if (! p1->tip) {
  2800.           tested += addTraverse(tr, p, p1->next->back,
  2801.                                 mintrav, maxtrav, bt, FALSE);
  2802.           tested += addTraverse(tr, p, p1->next->next->back,
  2803.                                 mintrav, maxtrav, bt, FALSE);
  2804.           }
  2805.  
  2806.         if (! p2->tip) {
  2807.           tested += addTraverse(tr, p, p2->next->back,
  2808.                                 mintrav, maxtrav, bt, FALSE);
  2809.           tested += addTraverse(tr, p, p2->next->next->back,
  2810.                                 mintrav, maxtrav, bt, FALSE);
  2811.           }
  2812.  
  2813.         hookup(p->next,       p1, p1z);  /*  Restore original tree */
  2814.         hookup(p->next->next, p2, p2z);
  2815.         initrav(tr->start);
  2816.         initrav(tr->start->back);
  2817.         }
  2818.       }   /* if (! p->tip) */
  2819.  
  2820. /* Moving subtree backward in tree.  Minimum move is 2 to avoid duplicates */
  2821.  
  2822.     q = p->back;
  2823.     if (! q->tip && maxtrav > 1) {
  2824.       q1 = q->next->back;
  2825.       q2 = q->next->next->back;
  2826.       if (! q1->tip && (!q1->next->back->tip || !q1->next->next->back->tip) ||
  2827.           ! q2->tip && (!q2->next->back->tip || !q2->next->next->back->tip)) {
  2828.         q1z = q1->z;
  2829.         q2z = q2->z;
  2830.         (void) removeNode(tr, q);
  2831.         cacheZ(tr);
  2832.         mintrav2 = mintrav > 2 ? mintrav : 2;
  2833.  
  2834.         if (! q1->tip) {
  2835.           tested += addTraverse(tr, q, q1->next->back,
  2836.                                 mintrav2 , maxtrav, bt, FALSE);
  2837.           tested += addTraverse(tr, q, q1->next->next->back,
  2838.                                 mintrav2 , maxtrav, bt, FALSE);
  2839.           }
  2840.  
  2841.         if (! q2->tip) {
  2842.           tested += addTraverse(tr, q, q2->next->back,
  2843.                                 mintrav2 , maxtrav, bt, FALSE);
  2844.           tested += addTraverse(tr, q, q2->next->next->back,
  2845.                                 mintrav2 , maxtrav, bt, FALSE);
  2846.           }
  2847.  
  2848.         hookup(q->next,       q1, q1z);  /*  Restore original tree */
  2849.         hookup(q->next->next, q2, q2z);
  2850.         initrav(tr->start);
  2851.         initrav(tr->start->back);
  2852.         }
  2853.       }   /* if (! q->tip && maxtrav > 1) */
  2854.  
  2855. /* Move other subtrees */
  2856.  
  2857.     if (! p->tip) {
  2858.       tested += rearrange(tr, p->next->back,       mintrav, maxtrav, bt);
  2859.       tested += rearrange(tr, p->next->next->back, mintrav, maxtrav, bt);
  2860.       }
  2861.  
  2862.     return  tested;
  2863.   } /* rearrange */
  2864.  
  2865. #ifdef SmallMemorySystem 
  2866. /* macintosh */
  2867. /* this is NOT define we want, but what Mac define fits all compilers !? */
  2868. long getpid(void) 
  2869. {
  2870.     return 0;
  2871. }
  2872. #endif
  2873.  
  2874. FILE *fopen_pid (filenm, mode, name_pid)
  2875.     char *filenm, *mode, *name_pid;
  2876.   { /* fopen_pid */
  2877.  
  2878.     (void) sprintf(name_pid, "%s.%d", filenm, getpid());
  2879.     return  fopen(name_pid, mode);
  2880.   } /* fopen_pid */
  2881.  
  2882.  
  2883. #if DeleteCheckpointFile
  2884. void  unlink_pid (filenm)
  2885.     char *filenm;
  2886.   { /* unlink_pid */
  2887.     char scr[512];
  2888.  
  2889.     (void) sprintf(scr, "%s.%d", filenm, getpid());
  2890.     unlink(scr);
  2891.   } /* unlink_pid */
  2892. #endif
  2893.  
  2894.  
  2895. void  writeCheckpoint (tr)
  2896.     tree  *tr;
  2897.   { /* writeCheckpoint */
  2898.     char   filename[128];
  2899.     FILE  *checkpointf;
  2900.     void   treeOut();
  2901.  
  2902.     checkpointf = fopen_pid(checkpointname, "a", filename);
  2903.     if (checkpointf) {
  2904.       treeOut(checkpointf, tr, 1);  /* 1 is for Newick format */
  2905.       (void) fclose(checkpointf);
  2906.       }
  2907.   } /* writeCheckpoint */
  2908.  
  2909.  
  2910. node * findAnyTip(p)
  2911.     nodeptr  p;
  2912.   { /* findAnyTip */
  2913.     return  p->tip ? p : findAnyTip(p->next->back);
  2914.   } /* findAnyTip */
  2915.  
  2916.  
  2917. void  optimize (tr, maxtrav, bt)
  2918.     tree     *tr;
  2919.     int       maxtrav;
  2920.     bestlist *bt;
  2921.   { /* optimize */
  2922.     nodeptr  p;
  2923.     int    mintrav, tested;
  2924.  
  2925.     if (tr->ntips < 4)  return;
  2926.  
  2927.     writeCheckpoint(tr);                    /* checkpoint the starting tree */
  2928.  
  2929.     if (maxtrav > tr->ntips - 3)  maxtrav = tr->ntips - 3;
  2930.     if (maxtrav <= tr->opt_level)  return;
  2931.  
  2932.     printf("      Doing %s rearrangements\n",
  2933.              (maxtrav == 1)            ? "local" :
  2934.              (maxtrav < tr->ntips - 3) ? "regional" : "global");
  2935.  
  2936.     /* loop while tree gets better  */
  2937.  
  2938.     do {
  2939.       (void) startOpt(bt, tr);
  2940.       mintrav = tr->opt_level + 1;
  2941.  
  2942.       /* rearrange must start from a tip or it will miss some trees */
  2943.  
  2944.       p = findAnyTip(tr->start);
  2945.       tested = rearrange(tr, p->back, mintrav, maxtrav, bt);
  2946.  
  2947. #     if Master
  2948.         getReturnedTrees(tr, bt, tested);
  2949. #     endif
  2950.  
  2951.       if (anerror)  return;
  2952.       bt->numtrees += tested;
  2953.       (void) setOptLevel(bt, maxtrav);
  2954.       (void) recallBestTree(bt, 1, tr);     /* recover best found tree */
  2955.  
  2956.       printf("      Tested %d alternative trees\n", tested);
  2957.       if (bt->improved) {
  2958.         printf("      Ln Likelihood =%14.5f\n", tr->likelihood);
  2959.         }
  2960.  
  2961.       writeCheckpoint(tr);                  /* checkpoint the new tree */
  2962.       } while (maxtrav > tr->opt_level);
  2963.  
  2964.   } /* optimize */
  2965.  
  2966.  
  2967. void coordinates (tr, p, lengthsum, tdptr)
  2968.     tree     *tr;
  2969.     nodeptr   p;
  2970.     double    lengthsum;
  2971.     drawdata *tdptr;
  2972.   { /* coordinates */
  2973.     /* establishes coordinates of nodes */
  2974.     double  x, z;
  2975.     nodeptr  q, first, last;
  2976.  
  2977.     if (p->tip) {
  2978.       p->xcoord = nint(over * lengthsum);
  2979.       p->ymax = p->ymin = p->ycoord = tdptr->tipy;
  2980.       tdptr->tipy += down;
  2981.       if (lengthsum > tdptr->tipmax) tdptr->tipmax = lengthsum;
  2982.       }
  2983.  
  2984.     else {
  2985.       q = p->next;
  2986.       do {
  2987.         z = q->z;
  2988.         if (z < zmin) z = zmin;
  2989.         x = lengthsum - fracchange * log(z);
  2990.         coordinates(tr, q->back, x, tdptr);
  2991.         q = q->next;
  2992.         } while (p == tr->start->back ? q != p->next : q != p);
  2993.  
  2994.       first = p->next->back;
  2995.       q = p;
  2996.       while (q->next != p) q = q->next;
  2997.       last = q->back;
  2998.       p->xcoord = nint(over * lengthsum);
  2999.       p->ycoord = (first->ycoord + last->ycoord)/2;
  3000.       p->ymin = first->ymin;
  3001.       p->ymax = last->ymax;
  3002.       }
  3003.   } /* coordinates */
  3004.  
  3005.  
  3006. void copyTrimmedName (cp1, cp2)
  3007.     char  *cp1, *cp2;
  3008.  { /* copyTrimmedName */
  3009.    char *ep;
  3010.  
  3011.    ep = cp1;
  3012.    while (*ep)  ep++;                              /* move forward to end */
  3013.    ep--;                                           /* move back to last */
  3014.    while (ep >= cp1 && white((int) *(ep)))  ep--;  /* trim white */
  3015.    while (cp1 <= ep)  *cp2++ = *cp1++;             /* copy to new end */
  3016.    *cp2 = 0;
  3017.  } /* copyTrimmedName */
  3018.  
  3019.  
  3020. void drawline (tr, i, scale)
  3021.     tree   *tr;
  3022.     int     i;
  3023.     double  scale;
  3024.     /* draws one row of the tree diagram by moving up tree */
  3025.     /* Modified to handle 1000 taxa, October 16, 1991 */
  3026.   { /* drawline */
  3027.     nodeptr  p, q, r, first, last;
  3028.     char    *nameptr;
  3029.     int  n, j, k, l, extra;
  3030.     boolean  done;
  3031.  
  3032.     p = q = tr->start->back;
  3033.     extra = 0;
  3034.  
  3035.     if (i == p->ycoord) {
  3036.       k = q->number - tr->mxtips;
  3037.       for (j = k; j < 1000; j *= 10) putchar('-');
  3038.       printf("%d", k);
  3039.       extra = 1;
  3040.       }
  3041.     else printf("   ");
  3042.  
  3043.     do {
  3044.       if (! p->tip) {
  3045.         r = p->next;
  3046.         done = FALSE;
  3047.         do {
  3048.           if ((i >= r->back->ymin) && (i <= r->back->ymax)) {
  3049.             q = r->back;
  3050.             done = TRUE;
  3051.             }
  3052.           r = r->next;
  3053.           } while (! done && (p == tr->start->back ? r != p->next : r != p));
  3054.  
  3055.         first = p->next->back;
  3056.         r = p;
  3057.         while (r->next != p) r = r->next;
  3058.         last = r->back;
  3059.         if (p == tr->start->back) last = p->back;
  3060.         }
  3061.  
  3062.       done = (p->tip) || (p == q);
  3063.       n = nint(scale*(q->xcoord - p->xcoord));
  3064.       if ((n < 3) && (! q->tip)) n = 3;
  3065.       n -= extra;
  3066.       extra = 0;
  3067.  
  3068.       if ((q->ycoord == i) && (! done)) {
  3069.         if (p->ycoord != q->ycoord) putchar('+');
  3070.         else                        putchar('-');
  3071.  
  3072.         if (! q->tip) {
  3073.           k = q->number - tr->mxtips;
  3074.           l = n - 3;
  3075.           for (j = k; j < 100; j *= 10)  l++;
  3076.           for (j = 1; j <= l; j++) putchar('-');
  3077.           printf("%d", k);
  3078.           extra = 1;
  3079.           }
  3080.         else for (j = 1; j <= n-1; j++) putchar('-');
  3081.         }
  3082.  
  3083.       else if (! p->tip) {
  3084.         if ((last->ycoord > i) && (first->ycoord < i) && (i != p->ycoord)) {
  3085.           putchar('!');
  3086.           for (j = 1; j <= n-1; j++) putchar(' ');
  3087.           }
  3088.         else for (j = 1; j <= n; j++) putchar(' ');
  3089.         }
  3090.  
  3091.       else
  3092.         for (j = 1; j <= n; j++) putchar(' ');
  3093.  
  3094.       p = q;
  3095.       } while (! done);
  3096.  
  3097.     if ((p->ycoord == i) && p->tip) {
  3098.       char  trimmed[nmlngth + 1];
  3099.  
  3100.       copyTrimmedName(p->name, trimmed);
  3101.       printf(" %s", trimmed);
  3102.       }
  3103.  
  3104.     putchar('\n');
  3105.   } /* drawline */
  3106.  
  3107.  
  3108. void printTree (tr)
  3109.     tree  *tr;
  3110.     /* prints out diagram of the tree */
  3111.   { /* printTree */
  3112.     drawdata  tipdata;
  3113.     double  scale;
  3114.     int  i, imax;
  3115.  
  3116.     if (treeprint) {
  3117.       putchar('\n');
  3118.       tipdata.tipy = 1;
  3119.       tipdata.tipmax = 0.0;
  3120.       coordinates(tr, tr->start->back, (double) 0.0, & tipdata);
  3121.       scale = 1.0 / tipdata.tipmax;
  3122.       imax = tipdata.tipy - down;
  3123.       for (i = 1; i <= imax; i++)  drawline(tr, i, scale);
  3124.       printf("\nRemember: ");
  3125.       if (outgropt) printf("(although rooted by outgroup) ");
  3126.       printf("this is an unrooted tree!\n\n");
  3127.       }
  3128.   } /* printTree */
  3129.  
  3130.  
  3131. double sigma (p, sumlrptr)
  3132.     nodeptr  p;
  3133.     double  *sumlrptr;
  3134.     /* compute standard deviation */
  3135.   { /* sigma */
  3136.     double  slope, sum, sumlr, z, zv, zz, lz,
  3137.             rat, suma, sumb, sumc, d2, d, li, temp, abzz, bczv, t3,
  3138.             fx1a, fx1c, fx1g, fx1t, fx1r, fx1y, fx2r, fx2y, w;
  3139.     double  *rptr;
  3140.     xtype   *x1a, *x1c, *x1g, *x1t, *x2a, *x2c, *x2g, *x2t;
  3141.     nodeptr  q;
  3142.     int  i, *wptr;
  3143.  
  3144.     q = p->back;
  3145.     while ((! p->x) || (! q->x)) {
  3146.       if (! (p->x)) newview(p);
  3147.       if (! (q->x)) newview(q);
  3148.       }
  3149.  
  3150.     x1a = &(p->x->a[0]);
  3151.     x1c = &(p->x->c[0]);
  3152.     x1g = &(p->x->g[0]);
  3153.     x1t = &(p->x->t[0]);
  3154.  
  3155.     x2a = &(q->x->a[0]);
  3156.     x2c = &(q->x->c[0]);
  3157.     x2g = &(q->x->g[0]);
  3158.     x2t = &(q->x->t[0]);
  3159.  
  3160.     z = p->z;
  3161.     if (z < zmin) z = zmin;
  3162.     lz = log(z);
  3163.  
  3164.     wptr = &(aliasweight[0]);
  3165.     rptr = &(ratvalue[0]);
  3166.     sum = sumlr = slope = 0.0;
  3167.  
  3168. #   ifdef Vectorize
  3169. #     pragma IVDEP
  3170. #   endif
  3171.  
  3172.     for (i = 0; i < endsite; i++) {
  3173.       rat  = *rptr++;
  3174.       zz   = exp(rat    * lz);
  3175.       zv   = exp(rat*xv * lz);
  3176.  
  3177.       fx1a = freqa * *x1a++;
  3178.       fx1g = freqg * *x1g++;
  3179.       fx1c = freqc * *x1c++;
  3180.       fx1t = freqt * *x1t++;
  3181.       fx1r = fx1a + fx1g;
  3182.       fx1y = fx1c + fx1t;
  3183.       suma = fx1a * *x2a + fx1c * *x2c + fx1g * *x2g + fx1t * *x2t;
  3184.       fx2r = freqa * *x2a++ + freqg * *x2g++;
  3185.       fx2y = freqc * *x2c++ + freqt * *x2t++;
  3186.       sumc = (fx1r + fx1y) * (fx2r + fx2y);
  3187.       sumb = fx1r * fx2r * invfreqr + fx1y * fx2y * invfreqy;
  3188.       abzz = zz * (suma - sumb);
  3189.       bczv = zv * (sumb - sumc);
  3190.       li = sumc + abzz + bczv;
  3191.       t3 = xv * bczv;
  3192.       d  = abzz + t3;
  3193.       d2 = rat * (abzz*(rat-1.0) + t3*(rat*xv-1.0));
  3194.  
  3195.       temp = rat * d / li;
  3196.       w = *wptr++;
  3197.       slope += w *  temp;
  3198.       sum   += w * (temp * temp - d2/li);
  3199.       sumlr += w * log(li/(suma+1.0E-300));
  3200.       }
  3201.  
  3202.     *sumlrptr = sumlr;
  3203.     return (sum > 1.0E-300) ? z*(-slope + sqrt(slope*slope + 3.841*sum))/sum
  3204.                             : 1.0;
  3205.   } /* sigma */
  3206.  
  3207.  
  3208. void describe (tr, p)
  3209.     tree    *tr;
  3210.     nodeptr  p;
  3211.     /* print out information for one branch */
  3212.   { /* describe */
  3213.     double  z, s, sumlr;
  3214.     nodeptr  q;
  3215.  
  3216.     q = p->back;
  3217.     printf("%4d          ", q->number - tr->mxtips);
  3218.     if (p->tip) printf("%s", p->name);
  3219.     else        printf("%4d      ", p->number - tr->mxtips);
  3220.  
  3221.     z = q->z;
  3222.     if (z <= zmin) printf("    infinity");
  3223.     else printf("%15.5f", -log(z)*fracchange);
  3224.  
  3225.     s = sigma(q, &sumlr);
  3226.     printf("     (");
  3227.     if (z + s >= zmax) printf("     zero");
  3228.     else printf("%9.5f", (double) -log(z + s)*fracchange);
  3229.     putchar(',');
  3230.     if (z - s <= zmin) printf("    infinity");
  3231.     else printf("%12.5f", (double) -log(z - s)*fracchange);
  3232.     putchar(')');
  3233.  
  3234.     if      (sumlr > 2.995 ) printf(" **");
  3235.     else if (sumlr > 1.9205) printf(" *");
  3236.     putchar('\n');
  3237.  
  3238.     if (! p->tip) {
  3239.       describe(tr, p->next->back);
  3240.       describe(tr, p->next->next->back);
  3241.       }
  3242.   } /* describe */
  3243.  
  3244.  
  3245. void summarize (tr)
  3246.     tree  *tr;
  3247.     /* print out branch length information and node numbers */
  3248.   { /* summarize */
  3249.     printf("Ln Likelihood =%14.5f\n", tr->likelihood);
  3250.     putchar('\n');
  3251.     printf(" Between        And             Length");
  3252.     printf("      Approx. Confidence Limits\n");
  3253.     printf(" -------        ---             ------");
  3254.     printf("      ------- ---------- ------\n");
  3255.  
  3256.     describe(tr, tr->start->back->next->back);
  3257.     describe(tr, tr->start->back->next->next->back);
  3258.     describe(tr, tr->start);
  3259.     putchar('\n');
  3260.     printf("     *  = significantly positive, P < 0.05\n");
  3261.     printf("     ** = significantly positive, P < 0.01\n\n\n");
  3262.   } /* summarize */
  3263.  
  3264.  
  3265. /*===========  This is a problem if tr->start->back is a tip!  ===========*/
  3266. /*  All routines should be contrived so that tr->start->back is not a tip */
  3267.  
  3268. char *treeString (treestr, tr, p, form)
  3269.     char  *treestr;
  3270.     tree  *tr;
  3271.     nodeptr  p;
  3272.     int  form;
  3273.     /* write string with representation of tree */
  3274.     /* form == 1 -> Newick tree */
  3275.     /* form == 2 -> Prolog fact */
  3276.   { /* treeString */
  3277.     double  x, z;
  3278.     char  *nameptr;
  3279.     int  n, c;
  3280.  
  3281.     if (p == tr->start->back) {
  3282.       if (form == 2) {
  3283.         (void) sprintf(treestr, "phylip_tree(");
  3284.         while (*treestr) treestr++;            /* move pointer to null */
  3285.         }
  3286.  
  3287.       (void) sprintf(treestr, "[&&%s: version = '%s'",
  3288.                                programName, programVersion);
  3289.       while (*treestr) treestr++;
  3290.  
  3291.       (void) sprintf(treestr, ", %s = %15.13g",
  3292.                                likelihood_key, tr->likelihood);
  3293.       while (*treestr) treestr++;
  3294.  
  3295.       (void) sprintf(treestr, ", %s = %d", ntaxa_key, tr->ntips);
  3296.       while (*treestr) treestr++;
  3297.  
  3298.       (void) sprintf(treestr,", %s = %d", opt_level_key, tr->opt_level);
  3299.       while (*treestr) treestr++;
  3300.  
  3301.       (void) sprintf(treestr, ", %s = %d", smoothed_key, tr->smoothed);
  3302.       while (*treestr) treestr++;
  3303.  
  3304.       (void) sprintf(treestr, "]%s", form == 2 ? ", " : " ");
  3305.       while (*treestr) treestr++;
  3306.       }
  3307.  
  3308.     if (p->tip) {
  3309.       *treestr++ = '\'';
  3310.       n = nmlngth;
  3311.       nameptr = p->name + nmlngth - 1;
  3312.       while (*nameptr-- == ' ' && n) n--;    /*  Trim trailing spaces */
  3313.       nameptr = p->name;
  3314.       while (n--) {
  3315.         if ((c = *nameptr++) == '\'')  *treestr++ = '\'';
  3316.         *treestr++ = c;
  3317.         }
  3318.       *treestr++ = '\'';
  3319.       }
  3320.  
  3321.     else {
  3322.       *treestr++ = '(';
  3323.       treestr = treeString(treestr, tr, p->next->back, form);
  3324.       *treestr++ = ',';
  3325.       treestr = treeString(treestr, tr, p->next->next->back, form);
  3326.       if (p == tr->start->back) {
  3327.         *treestr++ = ',';
  3328.         treestr = treeString(treestr, tr, p->back, form);
  3329.         }
  3330.       *treestr++ = ')';
  3331.       }
  3332.  
  3333.     if (p == tr->start->back) {
  3334.       (void) sprintf(treestr, ":0.0%s\n", (form != 2) ? ";" : ").");
  3335.       }
  3336.     else {
  3337.       z = p->z;
  3338.       if (z < zmin) z = zmin;
  3339.       x = -log(z) * fracchange;
  3340.       (void) sprintf(treestr, ": %8.6f", x);  /* prolog needs the space */
  3341.       }
  3342.  
  3343.     while (*treestr) treestr++;     /* move pointer up to null termination */
  3344.     return  treestr;
  3345.   } /* treeString */
  3346.  
  3347.  
  3348. void treeOut (treefile, tr, form)
  3349.     FILE  *treefile;
  3350.     tree  *tr;
  3351.     int    form;
  3352.     /* write out file with representation of final tree */
  3353.   { /* treeOut */
  3354.     int    c;
  3355.     char  *cptr, treestr[maxsp*(nmlngth+24)];
  3356.  
  3357.     (void) treeString(treestr, tr, tr->start->back, form);
  3358.     cptr = treestr;
  3359.     while (c = *cptr++) putc(c, treefile);
  3360.   } /* treeOut */
  3361.  
  3362.  
  3363. /*=======================================================================*/
  3364. /*                         Read a tree from a file                       */
  3365. /*=======================================================================*/
  3366.  
  3367.  
  3368. int treeFinishCom ()
  3369.   { /* treeFinishCom */
  3370.     int      ch;
  3371.     boolean  inquote;
  3372.  
  3373.     inquote = FALSE;
  3374.     while ((ch = getc(INFILE)) != EOF && (inquote || ch != ']')) {
  3375.       if (ch == '[' && ! inquote) {             /* comment; find its end */
  3376.         if ((ch = treeFinishCom()) == EOF)  break;
  3377.         }
  3378.       else if (ch == '\'') inquote = ! inquote;  /* start or end of quote */
  3379.       }
  3380.  
  3381.     return  ch;
  3382.   } /* treeFinishCom */
  3383.  
  3384.  
  3385. int treeGetCh ()
  3386.     /* get next nonblank, noncomment character */
  3387.   { /* treeGetCh */
  3388.     int  ch;
  3389.  
  3390.     while ((ch = getc(INFILE)) != EOF) {
  3391.       if (white(ch)) ;
  3392.       else if (ch == '[') {                   /* comment; find its end */
  3393.         if ((ch = treeFinishCom()) == EOF)  break;
  3394.         }
  3395.       else  break;
  3396.       }
  3397.  
  3398.     return  ch;
  3399.   } /* treeGetCh */
  3400.  
  3401.  
  3402. void  treeFlushLabel ()
  3403.   { /* treeFlushLabel */
  3404.     int      ch;
  3405.     boolean  done, quoted;
  3406.  
  3407.     if ((ch = treeGetCh()) == EOF)  return;
  3408.     done = (ch == ':' || ch == ',' || ch == ')'  || ch == '[' || ch == ';');
  3409.     if (! done && (quoted = (ch == '\'')))  ch = getc(INFILE);
  3410.  
  3411.     while (! done) {
  3412.       if (quoted) {
  3413.         if ((ch = findch('\'')) == EOF)  return;      /* find close quote */
  3414.         ch = getc(INFILE);                            /* check next char */
  3415.         if (ch != '\'') done = TRUE;                  /* not doubled quote */
  3416.         }
  3417.       else if (ch == ':' || ch == ',' || ch == ')'  || ch == '['
  3418.                          || ch == ';' || ch == '\n' || ch == EOF) {
  3419.         done = TRUE;
  3420.         }
  3421.       if (! done)  done = ((ch = getc(INFILE)) == EOF);
  3422.       }
  3423.  
  3424.     if (ch != EOF)  (void) ungetc(ch, INFILE);
  3425.   } /* treeFlushLabel */
  3426.  
  3427.  
  3428. int  findTipName (tr, ch)
  3429.     tree  *tr;
  3430.     int    ch;
  3431.   { /* findTipName */
  3432.     nodeptr  q;
  3433.     char  *nameptr, str[nmlngth+1];
  3434.     int  i, n;
  3435.     boolean  found, quoted, done;
  3436.  
  3437.     if (quoted = (ch == '\''))  ch = getc(INFILE);
  3438.     done = FALSE;
  3439.     i = 0;
  3440.  
  3441.     do {
  3442.       if (quoted) {
  3443.         if (ch == '\'') {
  3444.           ch = getc(INFILE);
  3445.           if (ch != '\'') done = TRUE;
  3446.           }
  3447.         else if (ch == EOF)
  3448.           done = TRUE;
  3449.         else if (ch == '\n' || ch == '\t')
  3450.           ch = ' ';
  3451.         }
  3452.       else if (ch == ':' || ch == ','  || ch == ')'  || ch == '['
  3453.                          || ch == '\n' || ch == EOF)
  3454.         done = TRUE;
  3455.       else if (ch == '_' || ch == '\t')
  3456.         ch = ' ';
  3457.  
  3458.       if (! done) {
  3459.         if (i < nmlngth)  str[i++] = ch;
  3460.         ch = getc(INFILE);
  3461.         }
  3462.       } while (! done);
  3463.  
  3464.     if (ch == EOF) {
  3465.       printf("ERROR: End-of-file in tree species name\n");
  3466.       return  0;
  3467.       }
  3468.  
  3469.     (void) ungetc(ch, INFILE);
  3470.     while (i < nmlngth)  str[i++] = ' ';     /*  Pad name */
  3471.  
  3472.     n = 1;
  3473.     do {
  3474.       q = tr->nodep[n];
  3475.       if (! (q->back)) {          /*  Only consider unused tips */
  3476.         i = 0;
  3477.         nameptr = q->name;
  3478.         do {found = str[i] == *nameptr++;}  while (found && (++i < nmlngth));
  3479.         }
  3480.       else
  3481.         found = FALSE;
  3482.       } while ((! found) && (++n <= tr->mxtips));
  3483.  
  3484.     if (! found) {
  3485.       i = nmlngth;
  3486.       do {str[i] = '\0';} while (i-- && (str[i] <= ' '));
  3487.       printf("ERROR: Cannot find data for tree species: %s\n", str);
  3488.       }
  3489.  
  3490.     return  (found ? n : 0);
  3491.   } /* findTipName */
  3492.  
  3493.  
  3494. double processLength ()
  3495.   { /* processLength */
  3496.     double  branch;
  3497.     int     ch;
  3498.     char    string[41];
  3499.  
  3500.     ch = treeGetCh();                            /*  Skip comments */
  3501.     if (ch != EOF)  (void) ungetc(ch, INFILE);
  3502.  
  3503.     if (fscanf(INFILE, "%lf", &branch) != 1) {
  3504.       printf("ERROR: Problem reading branch length in processLength:\n");
  3505.       if (fscanf(INFILE, "%40s", string) == 1)  printf("%s\n", string);
  3506.       anerror = TRUE;
  3507.       branch = 0.0;
  3508.       }
  3509.  
  3510.     return  branch;
  3511.   } /* processLength */
  3512.  
  3513.  
  3514. void  treeFlushLen ()
  3515.   { /* treeFlushLen */
  3516.     int  ch;
  3517.  
  3518.     if ((ch = treeGetCh()) == ':')
  3519.       (void) processLength();
  3520.     else if (ch != EOF)
  3521.       (void) ungetc(ch, INFILE);
  3522.  
  3523.   } /* treeFlushLen */
  3524.  
  3525.  
  3526. void  treeNeedCh (c1, where)
  3527.     int    c1;
  3528.     char  *where;
  3529.   { /* treeNeedCh */
  3530.     int  c2, i;
  3531.  
  3532.     if ((c2 = treeGetCh()) == c1)  return;
  3533.  
  3534.     printf("ERROR: Missing '%c' %s tree; ", c1, where);
  3535.     if (c2 == EOF) 
  3536.       printf("End-of-File");
  3537.     else {
  3538.       putchar('\'');
  3539.       for (i = 24; i-- && (c2 != EOF); c2 = getc(INFILE))  putchar(c2);
  3540.       putchar('\'');
  3541.       }
  3542.     printf(" found instead\n");
  3543.     anerror = TRUE;
  3544.   } /* treeNeedCh */
  3545.  
  3546.  
  3547. void  addElementLen (tr, p)
  3548.     tree    *tr;
  3549.     nodeptr  p;
  3550.   { /* addElementLen */
  3551.     double   z, branch;
  3552.     nodeptr  q;
  3553.     int      n, ch;
  3554.  
  3555.     if ((ch = treeGetCh()) == '(') {     /*  A new internal node */
  3556.       n = (tr->nextnode)++;
  3557.       if (n > 2*(tr->mxtips) - 2) {
  3558.         if (tr->rooted || n > 2*(tr->mxtips) - 1) {
  3559.           printf("ERROR: Too many internal nodes.  Is tree rooted?\n");
  3560.           printf("       Deepest splitting should be a trifurcation.\n");
  3561.           anerror = TRUE;
  3562.           return;
  3563.           }
  3564.         else {
  3565.           tr->rooted = TRUE;
  3566.           }
  3567.         }
  3568.       q = tr->nodep[n];
  3569.       addElementLen(tr, q->next);        if (anerror)  return;
  3570.       treeNeedCh(',', "in");             if (anerror)  return;
  3571.       addElementLen(tr, q->next->next);  if (anerror)  return;
  3572.       treeNeedCh(')', "in");             if (anerror)  return;
  3573.       treeFlushLabel();                  if (anerror)  return;
  3574.       }
  3575.  
  3576.     else {                               /*  A new tip */
  3577.       n = findTipName(tr, ch);
  3578.       if (n <= 0) {anerror = TRUE; return; }
  3579.       q = tr->nodep[n];
  3580.       if (tr->start->number > n)  tr->start = q;
  3581.       (tr->ntips)++;
  3582.       }                                  /* End of tip processing */
  3583.  
  3584.     if (tr->userlen) {
  3585.       treeNeedCh(':', "in");             if (anerror)  return;
  3586.       branch = processLength();          if (anerror)  return;
  3587.       z = exp(-branch / fracchange);
  3588.       if (z > zmax)  z = zmax;
  3589.       hookup(p, q, z);
  3590.       }
  3591.     else {
  3592.       treeFlushLen();                    if (anerror)  return;
  3593.       hookup(p, q, defaultz);
  3594.       }
  3595.   } /* addElementLen */
  3596.  
  3597.  
  3598. int saveTreeCom (comstrp)
  3599.     char  **comstrp;
  3600.   { /* saveTreeCom */
  3601.     int  ch;
  3602.     boolean  inquote;
  3603.  
  3604.     inquote = FALSE;
  3605.     while ((ch = getc(INFILE)) != EOF && (inquote || ch != ']')) {
  3606.       *(*comstrp)++ = ch;                        /* save character  */
  3607.       if (ch == '[' && ! inquote) {              /* comment; find its end */
  3608.         if ((ch = saveTreeCom(comstrp)) == EOF)  break;
  3609.         *(*comstrp)++ = ch;                      /* add ] */
  3610.         }
  3611.       else if (ch == '\'') inquote = ! inquote;  /* start or end of quote */
  3612.       }
  3613.  
  3614.     return  ch;
  3615.   } /* saveTreeCom */
  3616.  
  3617.  
  3618. boolean processTreeCom(tr)
  3619.     tree   *tr;
  3620.   { /* processTreeCom */
  3621.     int   ch, text_started, functor_read, com_open;
  3622.  
  3623.     /*  Accept prefatory "phylip_tree(" or "pseudoNewick("  */
  3624.  
  3625.     functor_read = text_started = 0;
  3626.     fscanf(INFILE, " p%nhylip_tree(%n", &text_started, &functor_read);
  3627.     if (text_started && ! functor_read) {
  3628.       fscanf(INFILE, "seudoNewick(%n", &functor_read);
  3629.       if (! functor_read) {
  3630.         printf("Start of tree 'p...' not understood.\n");
  3631.         anerror = TRUE;
  3632.         return;
  3633.         }
  3634.       }
  3635.  
  3636.     com_open = 0;
  3637.     fscanf(INFILE, " [%n", &com_open);
  3638.  
  3639.     if (com_open) {                              /* comment; read it */
  3640.       char  com[1024], *com_end;
  3641.  
  3642.       com_end = com;
  3643.       if (saveTreeCom(&com_end) == EOF) {        /* omits enclosing []s */
  3644.         printf("Missing end of tree comment.\n");
  3645.         anerror = TRUE;
  3646.         return;
  3647.         }
  3648.  
  3649.       *com_end = 0;
  3650.       (void) readKeyValue(com, likelihood_key, "%lg",
  3651.                                (void *) &(tr->likelihood));
  3652.       (void) readKeyValue(com, opt_level_key,  "%d",
  3653.                                (void *) &(tr->opt_level));
  3654.       (void) readKeyValue(com, smoothed_key,   "%d",
  3655.                                (void *) &(tr->smoothed));
  3656.  
  3657.       if (functor_read)  fscanf(INFILE, " ,");   /* remove trailing comma */
  3658.       }
  3659.  
  3660.     return (functor_read > 0);
  3661.   } /* processTreeCom */
  3662.  
  3663.  
  3664. void uprootTree (tr, p)
  3665.     tree   *tr;
  3666.     nodeptr p;
  3667.   { /* uprootTree */
  3668.     nodeptr  q, r, s;
  3669.     int  n;
  3670.  
  3671.     if (p->tip || p->back) {
  3672.       printf("ERROR: Unable to uproot tree.\n");
  3673.       printf("       Inappropriate node marked for removal.\n");
  3674.       anerror = TRUE;
  3675.       return;
  3676.       }
  3677.  
  3678.     n = --(tr->nextnode);               /* last internal node added */
  3679.     if (n != tr->mxtips + tr->ntips - 1) {
  3680.       printf("ERROR: Unable to uproot tree.  Inconsistent\n");
  3681.       printf("       number of tips and nodes for rooted tree.\n");
  3682.       anerror = TRUE;
  3683.       return;
  3684.       }
  3685.  
  3686.     q = p->next->back;                  /* remove p from tree */
  3687.     r = p->next->next->back;
  3688.     hookup(q, r, tr->userlen ? (q->z * r->z) : defaultz);
  3689.  
  3690.     q = tr->nodep[n];
  3691.     r = q->next;
  3692.     s = q->next->next;
  3693.     if (tr->ntips > 2 && p != q && p != r && p != s) {
  3694.       hookup(p,             q->back, q->z);   /* move connections to p */
  3695.       hookup(p->next,       r->back, r->z);
  3696.       hookup(p->next->next, s->back, s->z);
  3697.       }
  3698.  
  3699.     q->back = r->back = s->back = (nodeptr) NULL;
  3700.     tr->rooted = FALSE;
  3701.   } /* uprootTree */
  3702.  
  3703.  
  3704. void treeReadLen (tr)
  3705.     tree  *tr;
  3706.   { /* treeReadLen */
  3707.     nodeptr  p;
  3708.     int  i, ch;
  3709.     boolean  is_fact, found;
  3710.  
  3711.     for (i = 1; i <= tr->mxtips; i++) tr->nodep[i]->back = (node *) NULL;
  3712.     tr->start       = tr->nodep[tr->mxtips];
  3713.     tr->ntips       = 0;
  3714.     tr->nextnode    = tr->mxtips + 1;
  3715.     tr->opt_level   = 0;
  3716.     tr->log_f_valid = 0;
  3717.     tr->smoothed    = FALSE;
  3718.     tr->rooted      = FALSE;
  3719.  
  3720.     is_fact = processTreeCom(tr);
  3721.  
  3722.     p = tr->nodep[(tr->nextnode)++];
  3723.     treeNeedCh('(', "at start of");                  if (anerror)  return;
  3724.     addElementLen(tr, p);                            if (anerror)  return;
  3725.     treeNeedCh(',', "in");                           if (anerror)  return;
  3726.     addElementLen(tr, p->next);                      if (anerror)  return;
  3727.     if (! tr->rooted) {
  3728.       if ((ch = treeGetCh()) == ',') {        /*  An unrooted format */
  3729.         addElementLen(tr, p->next->next);            if (anerror)  return;
  3730.         }
  3731.       else {                                  /*  A rooted format */
  3732.         p->next->next->back = (nodeptr) NULL;
  3733.         tr->rooted = TRUE;
  3734.         if (ch != EOF)  (void) ungetc(ch, INFILE);
  3735.         }
  3736.       }
  3737.     treeNeedCh(')', "in");                           if (anerror)  return;
  3738.     treeFlushLabel();                                if (anerror)  return;
  3739.     treeFlushLen();                                  if (anerror)  return;
  3740.     if (is_fact) {
  3741.       treeNeedCh(')', "at end of");                  if (anerror)  return;
  3742.       treeNeedCh('.', "at end of");                  if (anerror)  return;
  3743.       }
  3744.     else {
  3745.       treeNeedCh(';', "at end of");                  if (anerror)  return;
  3746.       }
  3747.  
  3748.     if (tr->rooted)  uprootTree(tr, p->next->next);  if (anerror)  return;
  3749.     tr->start = p->next->next->back;  /* This is start used by treeString */
  3750.  
  3751.     initrav(tr->start);
  3752.     initrav(tr->start->back);
  3753.   } /* treeReadLen */
  3754.  
  3755.  
  3756. /*=======================================================================*/
  3757. /*                        Read a tree from a string                      */
  3758. /*=======================================================================*/
  3759.  
  3760.  
  3761. #if Master || Slave
  3762. int str_treeFinishCom (treestrp)
  3763.     char **treestrp;  /*  tree string pointer */
  3764.   { /* str_treeFinishCom */
  3765.     int ch;
  3766.     boolean  inquote;
  3767.  
  3768.     inquote = FALSE;
  3769.     while ((ch = *(*treestrp)++) != NULL && (inquote || ch != ']')) {
  3770.       if      (ch == '[' && ! inquote) {         /* comment; find its end */
  3771.         if ((ch = str_treeFinishCom(treestrp)) == NULL)  break;
  3772.         }
  3773.       else if (ch == '\'') inquote = ! inquote;  /* start or end of quote */
  3774.       }
  3775.     return  ch;
  3776.   } /* str_treeFinishCom */
  3777.  
  3778.  
  3779. int str_treeGetCh (treestrp)
  3780.     char **treestrp;  /*  tree string pointer */
  3781.     /* get next nonblank, noncomment character */
  3782.   { /* str_treeGetCh */
  3783.     int  ch;
  3784.     
  3785.     while ((ch = *(*treestrp)++) != NULL) {
  3786.       if (white(ch)) ;
  3787.       else if (ch == '[') {                   /* comment; find its end */
  3788.         if ((ch = str_treeFinishCom(treestrp)) == NULL)  break;
  3789.         }
  3790.       else  break;
  3791.       }
  3792.  
  3793.     return  ch;
  3794.   } /* str_treeGetCh */
  3795.  
  3796.  
  3797. void  str_treeFlushLabel (treestrp)
  3798.     char **treestrp;  /*  tree string pointer */
  3799.   { /* str_treeFlushLabel */
  3800.     int      ch;
  3801.     boolean  done, quoted;
  3802.  
  3803.     if ((ch = str_treeGetCh(treestrp)) == NULL)  done = TRUE;
  3804.     else {
  3805.       done = (ch == ':' || ch == ',' || ch == ')'  || ch == '[' || ch == ';');
  3806.       if (! done && (quoted = (ch == '\'')))  ch = *(*treestrp)++;
  3807.       }
  3808.  
  3809.     while (! done) {
  3810.       if (quoted) {
  3811.         if ((ch = str_findch(treestrp, '\'')) == NULL)  done = TRUE;
  3812.         else {
  3813.           ch = *(*treestrp)++;                        /* check next char */
  3814.           if (ch != '\'') done = TRUE;                /* not doubled quote */
  3815.           }
  3816.         }
  3817.       else if (ch == ':' || ch == ',' || ch == ')'  || ch == '['
  3818.                          || ch == ';' || ch == '\n' || ch == NULL) {
  3819.         done = TRUE;
  3820.         }
  3821.       if (! done)  done = ((ch = *(*treestrp)++) == NULL);
  3822.       }
  3823.  
  3824.     (*treestrp)--;
  3825.   } /* str_treeFlushLabel */
  3826.  
  3827.  
  3828. int  str_findTipName (treestrp, tr, ch)
  3829.     char **treestrp;  /*  tree string pointer */
  3830.     tree  *tr;
  3831.     int    ch;
  3832.   { /* str_findTipName */
  3833.     nodeptr  q;
  3834.     char  *nameptr, str[nmlngth+1];
  3835.     int  i, n;
  3836.     boolean  found, quoted, done;
  3837.  
  3838.     i = 0;
  3839.     if (quoted = (ch == '\'')) ch = *(*treestrp)++;
  3840.     done = FALSE;
  3841.  
  3842.     do {
  3843.       if (quoted) {
  3844.         if (ch == '\'') {
  3845.           ch = *(*treestrp)++;
  3846.           if (ch != '\'') done = TRUE;
  3847.           }
  3848.         else if (ch == NULL)
  3849.           done = TRUE;
  3850.         else if (ch == '\n' || ch == '\t')
  3851.           ch = ' ';
  3852.         }
  3853.       else if (ch == ':' || ch == ','  || ch == ')'  || ch == '['
  3854.                          || ch == '\n' || ch == NULL)
  3855.         done = TRUE;
  3856.       else if (ch == '_' || ch == '\t')
  3857.         ch = ' ';
  3858.  
  3859.       if (! done) {
  3860.         if (i < nmlngth)  str[i++] = ch;
  3861.         ch = *(*treestrp)++;
  3862.         }
  3863.       } while (! done);
  3864.  
  3865.     (*treestrp)--;
  3866.     if (ch == NULL) {
  3867.       printf("ERROR: NULL in tree species name\n");
  3868.       return  0;
  3869.       }
  3870.  
  3871.     while (i < nmlngth)  str[i++] = ' ';     /*  Pad name */
  3872.  
  3873.     n = 1;
  3874.     do {
  3875.       q = tr->nodep[n];
  3876.       if (! (q->back)) {          /*  Only consider unused tips */
  3877.         i = 0;
  3878.         nameptr = q->name;
  3879.         do {found = str[i] == *nameptr++;}  while (found && (++i < nmlngth));
  3880.         }
  3881.       else
  3882.         found = FALSE;
  3883.       } while ((! found) && (++n <= tr->mxtips));
  3884.  
  3885.     if (! found) {
  3886.       i = nmlngth;
  3887.       do {str[i] = '\0';} while (i-- && (str[i] <= ' '));
  3888.       printf("ERROR: Cannot find data for tree species: %s\n", str);
  3889.       }
  3890.  
  3891.     return  (found ? n : 0);
  3892.   } /* str_findTipName */
  3893.  
  3894.  
  3895. double str_processLength (treestrp)
  3896.     char **treestrp;   /*  tree string ponter */
  3897.   { /* str_processLength */
  3898.     double  branch;
  3899.     int     used;
  3900.  
  3901.     (void) str_treeGetCh(treestrp);                /*  Skip comments */
  3902.     (*treestrp)--;
  3903.  
  3904.     if (sscanf(*treestrp, "%lf%n", &branch, &used) != 1) {
  3905.       printf("ERROR: Problem reading branch length in str_processLength:\n");
  3906.       printf("%40s\n", *treestrp);
  3907.       anerror = TRUE;
  3908.       branch = 0.0;
  3909.       }
  3910.     else {
  3911.       *treestrp += used;
  3912.       }
  3913.  
  3914.     return  branch;
  3915.   } /* str_processLength */
  3916.  
  3917.  
  3918. void  str_treeFlushLen (treestrp)
  3919.     char **treestrp;   /*  tree string ponter */
  3920.   { /* str_treeFlushLen */
  3921.     int  ch;
  3922.  
  3923.     if ((ch = str_treeGetCh(treestrp)) == ':')
  3924.       (void) str_processLength(treestrp);
  3925.     else
  3926.       (*treestrp)--;
  3927.  
  3928.   } /* str_treeFlushLen */
  3929.  
  3930.  
  3931. void  str_treeNeedCh (treestrp, c1, where)
  3932.     char **treestrp;   /*  tree string pointer */
  3933.     int    c1;
  3934.     char  *where;
  3935.   { /* str_treeNeedCh */
  3936.     int  c2, i;
  3937.  
  3938.     if ((c2 = str_treeGetCh(treestrp)) == c1)  return;
  3939.  
  3940.     printf("ERROR: Missing '%c' %s tree; ", c1, where);
  3941.     if (c2 == NULL) 
  3942.       printf("NULL");
  3943.     else {
  3944.       putchar('\'');
  3945.       for (i = 24; i-- && (c2 != NULL); c2 = *(*treestrp)++)  putchar(c2);
  3946.       putchar('\'');
  3947.       }
  3948.     printf(" found instead\n");
  3949.     anerror = TRUE;
  3950.   } /* str_treeNeedCh */
  3951.  
  3952.  
  3953. void  str_addElementLen (treestrp, tr, p)
  3954.     char **treestrp;  /*  tree string pointer */
  3955.     tree    *tr;
  3956.     nodeptr  p;
  3957.   { /* str_addElementLen */
  3958.     double   z, branch;
  3959.     nodeptr  q;
  3960.     int      n, ch;
  3961.  
  3962.     if ((ch = str_treeGetCh(treestrp)) == '(') { /*  A new internal node */
  3963.       n = (tr->nextnode)++;
  3964.       if (n > 2*(tr->mxtips) - 2) {
  3965.         if (tr->rooted || n > 2*(tr->mxtips) - 1) {
  3966.           printf("ERROR: too many internal nodes.  Is tree rooted?\n");
  3967.           printf("Deepest splitting should be a trifurcation.\n");
  3968.           anerror = TRUE;
  3969.           return;
  3970.           }
  3971.         else {
  3972.           tr->rooted = TRUE;
  3973.           }
  3974.         }
  3975.       q = tr->nodep[n];
  3976.       str_addElementLen(treestrp, tr, q->next);        if (anerror)  return;
  3977.       str_treeNeedCh(treestrp, ',', "in");             if (anerror)  return;
  3978.       str_addElementLen(treestrp, tr, q->next->next);  if (anerror)  return;
  3979.       str_treeNeedCh(treestrp, ')', "in");             if (anerror)  return;
  3980.       str_treeFlushLabel(treestrp);                    if (anerror)  return;
  3981.       }
  3982.  
  3983.     else {                           /*  A new tip */
  3984.       n = str_findTipName(treestrp, tr, ch);
  3985.       if (n <= 0) {anerror = TRUE; return; }
  3986.       q = tr->nodep[n];
  3987.       if (tr->start->number > n)  tr->start = q;
  3988.       (tr->ntips)++;
  3989.       }                              /* End of tip processing */
  3990.  
  3991.     /*  Master and Slave always use lengths */
  3992.  
  3993.     str_treeNeedCh(treestrp, ':', "in");               if (anerror)  return;
  3994.     branch = str_processLength(treestrp);              if (anerror)  return;
  3995.     z = exp(-branch / fracchange);
  3996.     if (z > zmax)  z = zmax;
  3997.     hookup(p, q, z);
  3998.  
  3999.     return;
  4000.   } /* str_addElementLen */
  4001.  
  4002.  
  4003. boolean str_processTreeCom(tr, treestrp)
  4004.     tree   *tr;
  4005.     char  **treestrp;
  4006.   { /* str_processTreeCom */
  4007.     char  *com, *com_end;
  4008.     int  text_started, functor_read, com_open;
  4009.  
  4010.     com = *treestrp;
  4011.  
  4012.     functor_read = text_started = 0;
  4013.     sscanf(com, " p%nhylip_tree(%n", &text_started, &functor_read);
  4014.     if (functor_read) {
  4015.       com += functor_read;
  4016.       }
  4017.     else if (text_started) {
  4018.       com += text_started;
  4019.       sscanf(com, "seudoNewick(%n", &functor_read);
  4020.       if (! functor_read) {
  4021.         printf("Start of tree 'p...' not understood.\n");
  4022.         anerror = TRUE;
  4023.         return;
  4024.         }
  4025.       else {
  4026.         com += functor_read;
  4027.         }
  4028.       }
  4029.  
  4030.     com_open = 0;
  4031.     sscanf(com, " [%n", &com_open);
  4032.     com += com_open;
  4033.  
  4034.     if (com_open) {                              /* comment; read it */
  4035.     if (!(com_end = strchr(com, ']'))) {
  4036.         printf("Missing end of tree comment.\n");
  4037.         anerror = TRUE;
  4038.         return;
  4039.         }
  4040.  
  4041.       *com_end = 0;
  4042.       (void) readKeyValue(com, likelihood_key, "%lg",
  4043.                                (void *) &(tr->likelihood));
  4044.       (void) readKeyValue(com, opt_level_key,  "%d",
  4045.                                (void *) &(tr->opt_level));
  4046.       (void) readKeyValue(com, smoothed_key,   "%d",
  4047.                                (void *) &(tr->smoothed));
  4048.       *com_end = ']';
  4049.       com_end++;
  4050.  
  4051.       if (functor_read) {                          /* remove trailing comma */
  4052.         text_started = 0;
  4053.         sscanf(com_end, " ,%n", &text_started);
  4054.         com_end += text_started;
  4055.         }
  4056.  
  4057.       *treestrp = com_end;
  4058.       }
  4059.  
  4060.     return (functor_read > 0);
  4061.   } /* str_processTreeCom */
  4062.  
  4063.  
  4064. void str_treeReadLen (treestr, tr)
  4065.     char  *treestr;
  4066.     tree  *tr;
  4067.     /* read string with representation of tree */
  4068.   { /* str_treeReadLen */
  4069.     nodeptr  p;
  4070.     int  i;
  4071.     boolean  is_fact, found;
  4072.  
  4073.     for (i = 1; i <= tr->mxtips; i++) tr->nodep[i]->back = (node *) NULL;
  4074.     tr->start       = tr->nodep[tr->mxtips];
  4075.     tr->ntips       = 0;
  4076.     tr->nextnode    = tr->mxtips + 1;
  4077.     tr->opt_level   = 0;
  4078.     tr->log_f_valid = 0;
  4079.     tr->smoothed    = Master;
  4080.     tr->rooted      = FALSE;
  4081.  
  4082.     is_fact = str_processTreeCom(tr, &treestr);
  4083.  
  4084.     p = tr->nodep[(tr->nextnode)++];
  4085.     str_treeNeedCh(&treestr, '(', "at start of");    if (anerror)  return;
  4086.     str_addElementLen(&treestr, tr, p);              if (anerror)  return;
  4087.     str_treeNeedCh(&treestr, ',', "in");             if (anerror)  return;
  4088.     str_addElementLen(&treestr, tr, p->next);        if (anerror)  return;
  4089.     if (! tr->rooted) {
  4090.       if (str_treeGetCh(&treestr) == ',') {        /*  An unrooted format */
  4091.         str_addElementLen(&treestr, tr, p->next->next);
  4092.         if (anerror)  return;
  4093.         }
  4094.       else {                                       /*  A rooted format */
  4095.         p->next->next->back = (nodeptr) NULL;
  4096.         tr->rooted = TRUE;
  4097.         treestr--;
  4098.         }
  4099.       }
  4100.     str_treeNeedCh(&treestr, ')', "in");             if (anerror)  return;
  4101.     str_treeFlushLabel(&treestr);                    if (anerror)  return;
  4102.     str_treeFlushLen(&treestr);                      if (anerror)  return;
  4103.     if (is_fact) {
  4104.       str_treeNeedCh(&treestr, ')', "at end of");    if (anerror)  return;
  4105.       str_treeNeedCh(&treestr, '.', "at end of");    if (anerror)  return;
  4106.       }
  4107.     else {
  4108.       str_treeNeedCh(&treestr, ';', "at end of");    if (anerror)  return;
  4109.       }
  4110.  
  4111.     if (tr->rooted)  uprootTree(tr, p->next->next);  if (anerror)  return;
  4112.     tr->start = p->next->next->back;  /* This is start used by treeString */
  4113.  
  4114.     initrav(tr->start);
  4115.     initrav(tr->start->back);
  4116.   } /* str_treeReadLen */
  4117. #endif
  4118.  
  4119.  
  4120. void treeEvaluate (tr, bt)         /* Evaluate a user tree */
  4121.     tree     *tr;
  4122.     bestlist *bt;
  4123.   { /* treeEvaluate */
  4124.  
  4125.     if (Slave || ! tr->userlen)  smoothTree(tr, 4 * smoothings);
  4126.  
  4127.     (void) evaluate(tr, tr->start);
  4128.     if (anerror)  return;
  4129.  
  4130. #   if ! Slave
  4131.       (void) saveBestTree(bt, tr);
  4132. #   endif
  4133.   } /* treeEvaluate */
  4134.  
  4135.  
  4136. #if Master || Slave
  4137. FILE *freopen_pid (filenm, mode, stream)
  4138.     char *filenm, *mode;
  4139.     FILE *stream;
  4140.   { /* freopen_pid */
  4141.     char scr[512];
  4142.  
  4143.     (void) sprintf(scr, "%s.%d", filenm, getpid());
  4144.     return  freopen(scr, mode, stream);
  4145.   } /* freopen_pid */
  4146. #endif
  4147.  
  4148.  
  4149. void  showBestTrees (bt, tr, treefile)
  4150.     bestlist *bt;
  4151.     tree     *tr;
  4152.     FILE     *treefile;
  4153.   { /* showBestTrees */
  4154.     int     rank;
  4155.  
  4156.     for (rank = 1; rank <= bt->nvalid; rank++) {
  4157.       if (rank > 1) {
  4158.         if (rank != recallBestTree(bt, rank, tr))  break;
  4159.         }
  4160.       (void) evaluate(tr, tr->start);
  4161.       if (anerror)  return;
  4162.       if (tr->outgrnode->back)  tr->start = tr->outgrnode;
  4163.       printTree(tr);
  4164.       summarize(tr);
  4165.       if (treefile)  treeOut(treefile, tr, trout);
  4166.       }
  4167.   } /* showBestTrees */
  4168.  
  4169.  
  4170. #ifdef SmallMemorySystem
  4171. double  log_f0[maxpatterns];      /* Save a copy of best log_f */
  4172. #endif
  4173.  
  4174. void cmpBestTrees (bt, tr)
  4175.     bestlist *bt;
  4176.     tree     *tr;
  4177.   { /* cmpBestTrees */
  4178.     double  sum, sum2, sd, temp, bestscore;
  4179. #ifndef SmallMemorySystem
  4180.     double  log_f0[maxpatterns];      /* Save a copy of best log_f */
  4181. #endif
  4182.     double *log_f_ptr, *log_f0_ptr;
  4183.     int     i, j, num, besttips;
  4184.  
  4185.     num = bt->nvalid;
  4186.     if ((num <= 1) || (weightsum <= 1)) return;
  4187.  
  4188.     printf("Tree      Ln L        Diff Ln L       Its S.D.");
  4189.     printf("   Significantly worse?\n\n");
  4190.  
  4191.     for (i = 1; i <= num; i++) {
  4192.       if (i != recallBestTree(bt, i, tr))  break;
  4193.       if (! (tr->log_f_valid))  (void) evaluate(tr, tr->start);
  4194.  
  4195.       printf("%3d%14.5f", i, tr->likelihood);
  4196.       if (i == 1) {
  4197.         printf("  <------ best\n");
  4198.         besttips = tr->ntips;
  4199.         bestscore = tr->likelihood;
  4200.         log_f0_ptr = log_f0;
  4201.         log_f_ptr  = tr->log_f;
  4202.         for (j = 0; j < endsite; j++)  *log_f0_ptr++ = *log_f_ptr++;
  4203.         }
  4204.       else if (tr->ntips != besttips)
  4205.         printf("  (different number of species)\n");
  4206.       else {
  4207.         sum = sum2 = 0.0;
  4208.         log_f0_ptr = log_f0;
  4209.         log_f_ptr  = tr->log_f;
  4210.         for (j = 0; j < endsite; j++) {
  4211.           temp  = *log_f0_ptr++ - *log_f_ptr++;
  4212.           sum  += aliasweight[j] * temp;
  4213.           sum2 += aliasweight[j] * temp * temp;
  4214.           }
  4215.         sd = sqrt( weightsum * (sum2 - sum*sum/weightsum) / (weightsum-1) );
  4216.         printf("%14.5f%14.4f", tr->likelihood - bestscore, sd);
  4217.         printf("           %s\n", (sum > 1.95996 * sd) ? "Yes" : " No");
  4218.         }
  4219.       }
  4220.  
  4221.     printf("\n\n");
  4222.   } /* cmpBestTrees */
  4223.  
  4224.  
  4225. void makeUserTree (tr, bt)
  4226.     tree     *tr;
  4227.     bestlist *bt;
  4228.   { /* makeUserTree */
  4229.     char   filename[128];
  4230.     FILE  *treefile;
  4231.     int    nusertrees, which;
  4232.  
  4233.     if (fscanf(INFILE, "%d", &nusertrees) != 1 || findch('\n') == EOF) {
  4234.       printf("ERROR: Problem reading number of user trees\n");
  4235.       anerror = TRUE;
  4236.       return;
  4237.       }
  4238.  
  4239.     printf("User-defined %s:\n\n", (nusertrees == 1) ? "tree" : "trees");
  4240.  
  4241.     treefile = trout ? fopen_pid("treefile", "w", filename) : (FILE *) NULL;
  4242.  
  4243.     for (which = 1; which <= nusertrees; which++) {
  4244.       treeReadLen(tr);                     if (anerror)  return;
  4245.       treeEvaluate(tr, bt);                if (anerror)  return;
  4246.       if (tr->global <= 0) {
  4247.         if (tr->outgrnode->back)  tr->start = tr->outgrnode;
  4248.         printTree(tr);
  4249.         summarize(tr);
  4250.         if (treefile)  treeOut(treefile, tr, trout);
  4251.         }
  4252.       else {
  4253.         printf("%6d:  Ln Likelihood =%14.5f\n", which, tr->likelihood);
  4254.         }
  4255.       }
  4256.  
  4257.     if (tr->global > 0) {
  4258.       putchar('\n');
  4259.       if (! recallBestTree(bt, 1, tr)) anerror = TRUE;
  4260.       if (anerror)  return;
  4261.       printf("      Ln Likelihood =%14.5f\n", tr->likelihood);
  4262.       optimize(tr, tr->global, bt);
  4263.       if (anerror)  return;
  4264.       if (tr->outgrnode->back)  tr->start = tr->outgrnode;
  4265.       printTree(tr);
  4266.       summarize(tr);
  4267.       if (treefile)  treeOut(treefile, tr, trout);
  4268.       }
  4269.  
  4270.     if (treefile) {
  4271.       (void) fclose(treefile);
  4272.       printf("Tree also written to %s\n", filename);
  4273.       }
  4274.  
  4275.     putchar('\n');
  4276.     if (anerror)  return;
  4277.  
  4278.     cmpBestTrees(bt, tr);
  4279.   } /* makeUserTree */
  4280.  
  4281.  
  4282. #if Slave
  4283. void slaveTreeEvaluate (tr, bt)
  4284.     tree     *tr;
  4285.     bestlist *bt;
  4286.   { /* slaveTreeEvaluate */
  4287.     boolean done;
  4288.  
  4289.     do {
  4290.        requestForWork();
  4291.        receiveTree(&comm_master, tr);
  4292.        done = tr->likelihood > 0.0;
  4293.        if (! done) {
  4294.          treeEvaluate(tr, bt);         if (anerror) return;
  4295.          sendTree(&comm_master, tr);   if (anerror) return;
  4296.          }
  4297.        } while (! done);
  4298.   } /* slaveTreeEvaluate */
  4299. #endif
  4300.  
  4301.  
  4302. double randum (seed)
  4303.     longer  seed;
  4304.     /* random number generator */
  4305.   { /* randum */
  4306.     int  i, j, k, sum, mult[4];
  4307.     longer  newseed;
  4308.     double  x;
  4309.  
  4310.     mult[0] = 13;
  4311.     mult[1] = 24;
  4312.     mult[2] = 22;
  4313.     mult[3] =  6;
  4314.     for (i = 0; i <= 5; i++) newseed[i] = 0;
  4315.     for (i = 0; i <= 5; i++) {
  4316.       sum = newseed[i];
  4317.       k = MIN(i,3);
  4318.       for (j = 0; j <= k; j++) sum += mult[j] * seed[i-j];
  4319.       newseed[i] = sum;
  4320.       for (j = i; j <= 4; j++) {
  4321.         newseed[j+1] += newseed[j] >> 6;
  4322.         newseed[j] &= 63;
  4323.         }
  4324.       }
  4325.     newseed[5] &= 3;
  4326.     x = 0.0;
  4327.     for (i = 0; i <= 5; i++) x = 0.015625 * x + (seed[i] = newseed[i]);
  4328.     return  0.25 * x;
  4329.   } /* randum */
  4330.  
  4331.  
  4332. void makeDenovoTree (tr, bt)
  4333.     tree     *tr;
  4334.     bestlist *bt;
  4335.   { /* makeDenovoTree */
  4336.     char   filename[128];
  4337.     FILE  *treefile;
  4338.     nodeptr  p;
  4339.     int  enterorder[maxsp+1];      /*  random entry order */
  4340.     int  i, j, k, nextsp, newsp, maxtrav, tested;
  4341.  
  4342.     if (restart) {
  4343.       printf("Restarting from tree with the following sequences:\n");
  4344.       tr->userlen = TRUE;
  4345.       treeReadLen(tr);                        if (anerror)  return;
  4346.       smoothTree(tr, smoothings);             if (anerror)  return;
  4347.       (void) evaluate(tr, tr->start);         if (anerror)  return;
  4348.       (void) saveBestTree(bt, tr);            if (anerror)  return;
  4349.  
  4350.       for (i = 1, j = tr->ntips; i <= tr->mxtips; i++) { /* find loose tips */
  4351.         if (! tr->nodep[i]->back) enterorder[++j] = i;
  4352.         else {
  4353.           char  trimmed[nmlngth + 1];
  4354.  
  4355.           copyTrimmedName(tr->nodep[i]->name, trimmed);
  4356.           printf("   %s\n", trimmed);
  4357. #         if Master
  4358.             if (i>3) REPORT_ADD_SPECS;
  4359. #         endif
  4360.           }
  4361.         }
  4362.       putchar('\n');
  4363.       }
  4364.  
  4365.     else {                                           /* start from scratch */
  4366.       tr->ntips = 0;
  4367.       for (i = 1; i <= tr->mxtips; i++) enterorder[i] = i;
  4368.       }
  4369.  
  4370.     if (jumble) for (i = tr->ntips + 1; i <= tr->mxtips; i++) {
  4371.       j = randum(seed)*(tr->mxtips - tr->ntips) + tr->ntips + 1;
  4372.       k = enterorder[j];
  4373.       enterorder[j] = enterorder[i];
  4374.       enterorder[i] = k;
  4375.       }
  4376.  
  4377.     bt->numtrees = 1;
  4378.     if (tr->ntips < tr->mxtips)  printf("Adding species:\n");
  4379.  
  4380.     if (tr->ntips == 0) {
  4381.       for (i = 1; i <= 3; i++) {
  4382.         char  trimmed[nmlngth + 1];
  4383.  
  4384.         copyTrimmedName(tr->nodep[enterorder[i]]->name, trimmed);
  4385.         printf("   %s\n", trimmed);
  4386.         }
  4387.       tr->nextnode = tr->mxtips + 1;
  4388.       buildSimpleTree(tr, enterorder[1], enterorder[2], enterorder[3]);
  4389.       }
  4390.  
  4391.     if (anerror)  return;
  4392.  
  4393.     while (tr->ntips < tr->mxtips || tr->opt_level < tr->global) {
  4394.       maxtrav = (tr->ntips == tr->mxtips) ? tr->global : tr->partswap;
  4395.       if (maxtrav > tr->ntips - 3)  maxtrav = tr->ntips - 3;
  4396.  
  4397.       if (tr->opt_level >= maxtrav) {
  4398.         char  trimmed[nmlngth + 1];
  4399.  
  4400.         nextsp = ++(tr->ntips);
  4401.         newsp = enterorder[nextsp];
  4402.         p = tr->nodep[newsp];
  4403.  
  4404.         copyTrimmedName(p->name, trimmed);
  4405.         printf("   %s\n", trimmed);
  4406.  
  4407. #       if Master
  4408.           if (nextsp % DNAML_STEP_TIME_COUNT == 1) {
  4409.             REPORT_STEP_TIME;
  4410.             }
  4411.           REPORT_ADD_SPECS;
  4412. #       endif
  4413.  
  4414.         (void) buildNewTip(tr, p);
  4415.  
  4416.         resetBestTree(bt);
  4417.         cacheZ(tr);
  4418.         tested = addTraverse(tr, p->back, findAnyTip(tr->start)->back,
  4419.                              1, tr->ntips - 2, bt, fastadd);
  4420.         bt->numtrees += tested;
  4421.         if (anerror)  return;
  4422.  
  4423. #       if Master
  4424.           getReturnedTrees(tr, bt, tested);
  4425. #       endif
  4426.  
  4427.         printf("      Tested %d alternative trees\n", tested);
  4428.  
  4429.         (void) recallBestTree(bt, 1, tr);
  4430.         if (! tr->smoothed) {
  4431.           smoothTree(tr, smoothings);         if (anerror)  return;
  4432.           (void) evaluate(tr, tr->start);     if (anerror)  return;
  4433.           (void) saveBestTree(bt, tr);
  4434.           }
  4435.  
  4436.         if (tr->ntips == 4)  tr->opt_level = 1;  /* All 4 taxon trees done */
  4437.         maxtrav = (tr->ntips == tr->mxtips) ? tr->global : tr->partswap;
  4438.         if (maxtrav > tr->ntips - 3)  maxtrav = tr->ntips - 3;
  4439.         }
  4440.  
  4441.       printf("      Ln Likelihood =%14.5f\n", tr->likelihood);
  4442.       optimize(tr, maxtrav, bt);
  4443.  
  4444.       if (anerror)  return;
  4445.       }
  4446.  
  4447.     printf("\nExamined %d %s\n", bt->numtrees,
  4448.                                  bt->numtrees != 1 ? "trees" : "tree");
  4449.  
  4450.     treefile = trout ? fopen_pid("treefile", "w", filename) : (FILE *) NULL;
  4451.     showBestTrees(bt, tr, treefile);
  4452.     if (treefile) {
  4453.       (void) fclose(treefile);
  4454.       printf("Tree also written to %s\n\n", filename);
  4455.       }
  4456.  
  4457.     cmpBestTrees(bt, tr);
  4458.  
  4459. #   if DeleteCheckpointFile
  4460.       unlink_pid(checkpointname);
  4461. #   endif
  4462.   } /* makeDenovoTree */
  4463.  
  4464. /*==========================================================================*/
  4465. /*                             "main" routine                               */
  4466. /*==========================================================================*/
  4467.  
  4468. #ifdef SmallMemorySystem
  4469. static tree      curtree;   /*  current tree */
  4470. static bestlist  bestree;   /*  topology of best found tree */
  4471. #endif
  4472.  
  4473. #ifdef MACINTOSH
  4474. /* for use as child application, dgg, march'95 */
  4475. int RealMain(int argc, char** argv)
  4476. #else
  4477. #if Sequential
  4478.   main ()
  4479. #else
  4480.   slave ()
  4481. #endif
  4482. #endif
  4483.  
  4484.   { /* DNA Maximum Likelihood */
  4485. #   if Master
  4486.       int starttime, inputtime, endtime;
  4487. #   endif
  4488.  
  4489. #   if Sequential
  4490.         /* dgg addition */
  4491.       long starttime, inputtime, endtime;
  4492.       double millisecs;
  4493. #   endif
  4494.  
  4495. #   if Master || Slave
  4496.       int my_id, nprocs, type, from, sz;
  4497.       char *msg;
  4498. #   endif
  4499.  
  4500. #ifndef SmallMemorySystem
  4501.     tree      curtree;   /*  current tree */
  4502.     bestlist  bestree;   /*  topology of best found tree */
  4503. #endif
  4504.     tree      *tr;   /*  current tree */
  4505.     bestlist  *bt;   /*  topology of best found tree */
  4506.  
  4507. #   if Debug
  4508.       {
  4509.         char debugfilename[128];
  4510.         debug = fopen_pid("dnaml_debug", "w", debugfilename);
  4511.         }
  4512. #   endif
  4513.  
  4514. #   if Master
  4515.       starttime = p4_clock();
  4516.       nprocs = p4_num_total_slaves();
  4517.  
  4518.       if ((OUTFILE = freopen_pid("master.out", "w", stdout)) == NULL) {
  4519.         fprintf(stderr, "Could not open output file\n");
  4520.         exit(1);
  4521.         }
  4522.  
  4523.       /* Receive input file name from host */
  4524.       type = DNAML_FILE_NAME;
  4525.       from = DNAML_HOST_ID;
  4526.       msg  = NULL;
  4527.       p4_recv(&type, &from, &msg, &sz);
  4528.       if ((INFILE = fopen(msg, "r")) == NULL) {
  4529.         fprintf(stderr, "master could not open input file %s\n", msg);
  4530.         exit(1);
  4531.         }
  4532.       p4_msg_free(msg);
  4533.  
  4534.       open_link(&comm_slave);
  4535. #   endif
  4536.  
  4537. #   if Sequential
  4538.       starttime = clock();
  4539. #ifdef NoSTDIN
  4540.             {
  4541.             char *fname = "infile";
  4542.       if ((INFILE = fopen(fname, "r")) == NULL) 
  4543.       do {
  4544.           char filename[512];
  4545.           filename[0]= 0;
  4546.         fprintf(stderr, "Could not open input file '%s'\n", fname);
  4547.         fprintf(stderr, "Enter path-name for data file: ");
  4548.         scanf( "%s", filename);
  4549.         fname= filename;
  4550.            INFILE = fopen(fname, "r");
  4551.       } while (!INFILE && *fname > ' ');
  4552.       if (INFILE == NULL) {
  4553.              fprintf(stderr, "Could not open input file '%s'\n", fname);
  4554.             exit(1);
  4555.           }
  4556.       }
  4557. #endif
  4558. #        endif
  4559.  
  4560. #  if DNAML_STEP
  4561.       begin_step_time = starttime;
  4562. #  endif
  4563.  
  4564. #   if Slave
  4565.       my_id = p4_get_my_id();
  4566.       nprocs = p4_num_total_slaves();
  4567.  
  4568.       /* Receive input file name from host */
  4569.       type = DNAML_FILE_NAME;
  4570.       from = DNAML_HOST_ID;
  4571.       msg  = NULL;
  4572.       p4_recv(&type, &from, &msg, &sz);
  4573.       if ((INFILE = fopen(msg, "r")) == NULL) {
  4574.         fprintf(stderr, "slave could not open input file %s\n",msg);
  4575.         exit(1);
  4576.         }
  4577.       p4_msg_free(msg);
  4578.  
  4579. #     ifdef P4DEBUG
  4580.         if ((OUTFILE = freopen_pid("slave.out", "w", stdout)) == NULL) {
  4581.           fprintf(stderr, "Could not open output file\n");
  4582.           exit(1);
  4583.           }
  4584. #     else
  4585.         if ((OUTFILE = freopen("/dev/null", "w", stdout)) == NULL) {
  4586.           fprintf(stderr, "Could not open output file\n");
  4587.           exit(1);
  4588.           }
  4589. #     endif
  4590.  
  4591.       open_link(&comm_master);
  4592. #   endif
  4593.  
  4594.     anerror = FALSE;
  4595.     tr = & curtree;
  4596.     bt = & bestree;
  4597.     bt->ninit = 0;
  4598.     getinput(tr); 
  4599.     if (anerror)  return 1;
  4600.  
  4601. #   if Master
  4602.       inputtime = p4_clock();
  4603.       printf("Input time %d milliseconds\n",inputtime-starttime);
  4604.       REPORT_STEP_TIME;
  4605. #   endif
  4606.  
  4607. #   if Slave
  4608.       (void) fclose(INFILE);
  4609. #   endif
  4610.  
  4611. #   if Sequential
  4612.       inputtime = clock();
  4613.       millisecs= TIMETOSECS(inputtime - starttime);
  4614.       printf("Input time %6.4f seconds\n",millisecs);
  4615. #   ifdef NoSTDIN
  4616.       (void) fclose(INFILE);
  4617. #   endif
  4618. #   endif
  4619.  
  4620. /*  The material below would be a loop over jumbles and/or boots */
  4621.  
  4622.     makeweights();                              if (anerror)  return 1;
  4623.     makevalues(tr);                             if (anerror)  return 1;
  4624.     if (freqsfrom) {
  4625.       empiricalfreqs(tr);                       if (anerror)  return 1;
  4626.       getbasefreqs();
  4627.       }
  4628.  
  4629.     linkxarray(3, 3, & freextip, & usedxtip);   if (anerror)  return 1;
  4630.     setupnodex(tr);                             if (anerror)  return 1;
  4631.  
  4632. #   if Slave
  4633.       slaveTreeEvaluate(tr, bt);
  4634. #   else
  4635.       (void) initBestTree(bt, nkeep, endsite);      if (anerror)  return 1;
  4636.       if (! tr->usertree) {makeDenovoTree(tr, bt);  if (anerror)  return 1;}
  4637.       else                {makeUserTree(tr, bt);    if (anerror)  return 1;}
  4638.       freeBestTree(bt);                             if (anerror)  return 1;
  4639. #   endif
  4640.  
  4641. /*  Endpoint for jumble and/or boot loop */
  4642.  
  4643. #   if Master
  4644.       tr->likelihood = 1.0;             /* terminate slaves */
  4645.       sendTree(&comm_slave, tr);
  4646. #   endif
  4647.  
  4648.     freeTree(tr);
  4649.  
  4650. #   if Master
  4651.       close_link(&comm_slave);
  4652.       (void) fclose(INFILE);
  4653.  
  4654.       REPORT_STEP_TIME;
  4655.       endtime = p4_clock();
  4656.       printf("Execution time %d milliseconds\n", endtime - inputtime);
  4657.       (void) fclose(OUTFILE);
  4658. #   endif
  4659.  
  4660. #   if Slave
  4661.       close_link(&comm_master);
  4662.       (void) fclose(OUTFILE);
  4663. #   endif
  4664.  
  4665. #   if Sequential
  4666.         /* dgg addition */
  4667.       endtime = clock();
  4668.       millisecs= TIMETOSECS(endtime - inputtime);
  4669.       printf("Execution time %6.4f seconds\n", millisecs);
  4670. #   endif
  4671.  
  4672. #   if Debug
  4673.       (void) fclose(debug);
  4674. #   endif
  4675.  
  4676. #   if Master || Slave
  4677.       p4_send(DNAML_DONE, DNAML_HOST_ID, NULL, 0);
  4678. #   else
  4679.       return 0;
  4680. #   endif
  4681.   } /* DNA Maximum Likelihood */
  4682.